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CHAPTER  1 

INTRODUCTION 

With  the  continuous  advances  in  solid  state  device  technology  in  terms  of  both 
speed  and  size,  a  number  of  challenges  are  presented  to  the  packaging  engineer  who  has  to 
accommodate  a  multitude  of  these  devices  on  the  limited  real  estate  available  on  the  circuit 
board.  In  the  absence  of  reliable  design  tools,  digital  circuit  packaging  is  carried  out 
through  a  costly  and  lengthy  trial  and  error  process.  Thus,  a  critical  need  exists  for 
computer-aided  design  tools  that  are  capable  of  predicting  the  electrical  performance  of  a 
given  package,  e.g.,  crosstalk  and  signal  distortion,  in  a  reliable  fashion.  Understanding 
how  printed  circuit  layouts  can  affect  these  quantities  is  essential  for  improved  package 
designs. 

Similarly,  in  the  area  of  aircraft  design  the  interaction  of  the  aircraft  with 
electromagnetic  waves  is  of  crucial  importance.  In  recent  years,  the  reduction  of  the 
amount  of  radar  power  reflected  by  a  complex  object  like  an  aircraft,  known  as  the  radar 
cross  section  (RCS),  has  become  a  major  objective  of  aircraft  designers.  Therefore,  there 
is  a  critical  need  for  efficient  electromagnetic  computer-aided  design  tools  in  the  area  of 
circuit  design,  where  the  major  concern  is  the  reduction  of  crosstalk  and  signal  distortion 
for  a  printed  circuit  board,  as  well  as  in  the  area  of  aircraft  design,  where  the  major  concern 
is  the  reduction  of  radar  cross  section  of  an  aircraft. 

A  plethora  of  techniques  is  available  in  the  literature  for  dealing  with  different 
transmission  line  and  electromagnetic  scattering  problems,  e.g.,  the  method  of  moments, 
Fourier  transform  method,  variational  method,  and  conformal  mapping  method.  However, 
all  of  these  approaches  are  limited  in  their  application  to  homogeneous  media  and  simple 
geometries.  Furthermore,  for  large  problems,  these  methods  become  numerically 
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inefficient  because  they  give  rise  to  full  matrices  which  require  large  computer  memory  and 
solution  times,  especially  for  open  region  problems.  In  contrast,  the  finite  element  method 
(FEM)  can  handle  complex-shaped  structures  and  highly  inhomogeneous  dielectrics.  In  the 
FEM,  the  region  of  interest  is  bounded  by  an  artificial  boundary  to  limit  the  number  of 
unknowns.  Over  the  bounded  region,  the  Helmholtz  or  Laplace  equation  is  solved  at  a 
finite  number  of  grid  points.  These  equations  are  discretized  through  the  use  of  a  weak 
variational  formulation.  Although  the  resulting  matrix  equation  is  substantially  larger  than 
that  obtained  in  other  methods,  especially  the  method  of  moments,  it  is,  nevertheless, 
highly  sparse  and  can  be  inverted  rather  efficiently  using  special  algorithms  that  store  only 
the  nonzero  entries  of  the  matrix.  The  major  advantage  in  using  FEM  is  the  simplicity  with 
which  complex-shaped  structures  can  be  modeled.  Another  advantage  is  the  great 
improvement  in  storage  and  even  in  the  computational  time,  over  the  method  of  moments, 
especially  in  the  case  of  inhomogeneous  dielectric  scatterers  where  a  volume  formulation  of 
the  integral  equation  is  required  These  features  of  the  FEM  formulation  make  it  a  good 
candidate  for  CAD  packages.  However,  one  drawback  of  FEM  is  that,  in  dealing  with 
open  region  problems,  they  require  the  introduction  of  an  artificial  outer  boundary  in  order 
to  limit  the  number  of  node  points  to  a  manageable  size.  The  major  difficulty  encountered 
when  using  FEM  is  how  to  find  the  proper  boundary  condition  that  can  be  applied  on  the 
artificial  outer  boundary  to  make  it  as  transparent  as  possible  in  other  words,  how  to 
impose  the  corresponding  behavior  at  infinity  on  the  finite  distance  boundary  and  obtain  an 
accurate  solution  in  the  interior  region.  The  answer  to  this  important  question  is  of  a  crucial 
importance  for  the  development  of  efficient  electromagnetic  computer-aided  design  tools. 
There  have  been  different  approaches  to  model  the  outer  boundary  both  for  electromagnetic 
scattering  and  digital  circuit  problems. 

For  the  quasi-static  analysis  of  digital  circuit  problems,  it  is  customary  to  use  either 
the  natural  boundary  condition  or  the  infinite  elements  for  mesh  truncation.  The  first 
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approach  consists  of  introducing  a  fictitious  conducting  enclosure  at  the  outer  boundary 
[l]-[4].  This  approach  gives  satisfactory  results  only  if  the  actual  field  has  decayed 
sufficiently  well  as  it  reaches  the  outer  boundary.  In  most  cases,  this  approach,  which 
assumes  that  the  field  decays  significandy  before  reaching  the  outer  boundary,  results  in  an 
undesirably  large  mesh,  especially  for  three-dimensional  problems.  The  second  approach 
uses  the  infinite  elements  which  extend  to  infinity  [5]-[7].  Although  this  approach  is 
superior  to  the  artificial  p.e.c.  boundary  method,  it  also  has  its  own  drawbacks.  First,  the 
infinite  elements  require  special  care  during  the  matrix  filling.  Second,  one  needs  to 
assume  a  certain  field  behavior  within  the  infinite  elements,  and  this  behavior  may  not  be 
known. 

For  electromagnetic  scattering  problems,  the  Helmholtz  type  of  equation  is  solved. 
In  order  to  model  the  field  behavior  at  the  outer  boundary,  the  so-called  local  and  local 
boundary  conditions  have  been  used.  In  the  first  type,  the  surface  integral  equation 
involving  the  free-space  Green’s  function  is  used  as  the  boundary  constraint  [8]-[l  1]  and 
[13]-[15].  Mei  [12]  has  used  a  similar  approach  where  an  eigenfunction  expansion  of  the 
field  and  its  derivative  are  used  as  constraints  on  the  outer  boundary.  Both  of  these 
methods  are  sometimes  classified  in  the  literature  under  the  umbrella  of  global  or  local 
boundary  conditions.  The  reason  for  this  classification  is  that  all  of  the  field  values  at  the 
boundary  are  related  through  the  surface  integral  constraint.  An  important  drawback  in 
using  a  nonlocal  type  of  boundary  condition  is  that  it  spoils  the  sparsity  of  the  system 
matrix  and,  consequently,  becomes  time-consuming  for  large  problems. 

In  the  local  boundary  condition  approach,  an  asymptotic  differential  boundary 
operator,  referred  to  in  the  literature  as  an  Absorbing  Boundary  Condition  (ABC)  [16]- 
[28],  is  used  as  a  constraint  on  the  normal  derivative  of  the  field  at  the  outer  boundary. 
This  boundary  operator  attempts  to  impose  an  outward  propagating  character  on  the 
scattered  field,  i.e.,  tries  to  eliminate  the  nonphysical  reflections  at  the  boundary  that 
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generate  incoming  waves.  The  ABC  operator  is  only  asymptotic  in  nature,  and  does  not 
satisfy  the  near-field  radiation  condition  in  the  exact  sense.  Thus,  it  introduces  an  error  in 
the  FEM  solution  because  it  is  not  totally  free  of  incoming  waves  generated  by  the  artificial 
outer  boundary;  however,  this  error  is  often  quite  small,  and  hence  acceptable  for  many 
practical  applications.  Furthermore,  in  contrast  to  the  local  boundary  conditions,  the  ABC 
does  preserve  the  sparsity  of  the  discretized  FEM  matrix  and  is,  therefore,  an  attractive 
candidate  for  numerical  applicadons.  Nevertheless,  there  have  been  two  major  limitations 
associated  with  this  type  of  boundary  conditions.  First,  most  attempts  to  solve  the  open 
region  scattering  problem  using  the  local-type  boundary  operators  have  been  limited  to 
separable  types  of  geometry  for  the  outer  boundary.  This  may  be  attributable  to  the  form  in 
which  these  operators  were  cast  in  their  original  versions.  However,  for  a  class  of 
scatterers  that  are  long,  such  as  a  perfectly  conducting  strip  or  an  air  foil,  truncating  the 
open  region  with  a  circular  outer  boundary  requires  solving  for  the  field  over  a  very  large 
mesh  region  surrounding  the  scatterer,  which,  in  turn,  requires  the  solution  of  a  large 
matrix.  In  such  situations,  the  FEM  becomes  less  efficient  than  the  method  of  moments. 
Second,  most  forms  of  the  absorbing  boundary  condition  operators  have  been  based  on  the 
use  of  only  the  first  few  terms  of  the  asymptotic  representation  of  the  solution  to  the 
differential  equation.  In  [24]  and  r27]-[28],  it  was  demonstrated  that  while  the  absorbing 
boundaiy  condition,  which  is  based  on  the  first  terms  of  the  series,  works  quite  well  for  the 
lower-order  harmonics,  it  exhibits  a  significant  error  for  the  higher-order  harmonics, 
especially  if  the  outer  boundary  is  placed  very  close  to  the  surface  of  the  object. 

The  goal  of  this  study  is  the  circumvent  some  of  the  above-mentioned  problems  that 
are  encountered  in  FEM  mesh  truncation.  First,  the  Bayliss,  Gunzburger,  and  Turkel 
(BGT)  [9]  boundary  condition  is  generalized  so  as  to  make  it  applicable  to  an  arbitrary, 
rather  than  circular,  outer  boundary.  The  use  of  the  generalized  version  of  the  BGT 
enables  one  to  reduce  the  number  of  node  points  significantly  and  solve  larger  sized 
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problems  than  had  been  possible  in  the  past.  Second,  an  asymptotic  differential  operator 
for  the  quasi-static  analysis  of  two-dimensional  transmission  line  structures  is  obtained. 
This  operator,  which  is  also  valid  for  an  arbitrary  outer  boundary,  not  only  allows  one  to 
bring  the  outer  boundary  much  closer  than  is  possible  with  a  p.e.c.  artificial  boundary,  but 
does  not  suffer  from  the  complications  associated  with  the  infinite  elements.  Third,  a  set  of 
asymptotic  boundary  condition  operators  for  finite  element  quasi-TEM  analysis  of  three- 
dimensional  transmission  line  discontinuities  is  derived  from  the  general  solution  to  the 
three-dimensional  Laplace  equation  in  spherical  coordinates.  The  second-order  three- 
dimensional  asymptotic  boundary  condition  is  then  applied  on  a  box-shaped  outer 
boundary  for  the  purpose  of  truncating  the  mesh  in  an  efficient  manner.  With  'his 
boundary  condition,  it  becomes  possible  to  consider  general  three-dimensional 
discontinuities  and  consider  more  practical  problems.  Fourth,  the  general  form  of  the 
solution  to  the  two-dimensional  Laplace  equation  is  used  to  derive  a  higher-order 
asymptotic  boundary  condition  for  transmission  line  circuits.  This  boundary  condition, 
unlike  the  one  discussed  earlier  which  assumes  that  in  the  far  region  the  solution  can 
adequately  be  represented  by  the  first  two  terms  of  the  series,  requires  that  the  asymptotic 
representation  be  a  combination  of  lower-  and  higher-order  terms.  Therefore,  it  corrects,  to 
a  good  degree,  the  eTor  caused  by  the  neglecting  of  the  higher-order  terms  by  the  simple 
asymptotic  boundary  condition,  which,  in  turn,  yields  a  significant  improvement  in  the 
finite  element  solution.  Last,  a  higher-order  absorbing  boundary  condition  for  scattering 
problems  is  derived  from  the  asymptotic  solution  of  the  Helmholtz  equation.  This  higher- 
order  boundary  condition,  which  takes  into  account  both  the  lower-  and  higher-order 
harmonics,  considerably  reduces  the  error  caused  by  the  neglecting  of  the  higher-order 
harmonics  by  most  of  the  available  absorbing  boundary  conditions. 

Chapter  2  gives  some  background  on  the  finite  element  method  and  the  absorbing 
and  asymptotic  boundary  conditions.  Chapter  3  deals  with  the  open  region  scattering 
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problem  and  the  generalization  of  the  original  BGT  boundary  condition  so  as  to  make  it 
applicable  to  an  arbitrary  outer  boundary.  The  asymptotic  boundary  condition  for  finite 
element  quasi-TEM  analysis  of  two-dimensional  transmission  line  structures  is  discussed  in 
Chapter  4.  The  corresponding  analysis  for  three-dimensional  transmission  line 
discontinuities  is  addressed  in  Chapter  5.  Chapter  6  deals  with  the  higher-order  boundary 
condition  for  both  the  scattering  as  well  as  the  transmission  line  circuit  problems  and  the 
improvement  on  the  simple  boundary  conditions.  In  Chapter  7  the  conclusions  drawn  from 
this  study  are  discussed. 
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CHAPTER  2 

THE  FINITE  ELEMENT  METHOD,  THE  ASYMPTOTIC, 

AND  THE  ABSORBING  BOUNDARY  CONDITIONS 

2.1  Introduction 

The  finite  element  method  (FEM)  has  become  an  important  numerical  technique  to 
solve  open  region  problems  because  of  its  flexibility  in  handling  arbitrary  structures  and 
highly  inhomogeneous  materials.  However,  it  must  deal  with  the  practical  problem  of 
mesh  truncation  and  the  large  number  of  mesh  nodes.  By  imposing  an  asymptotic  or  an 
absorbing  boundary  condition  (ABC)  on  the  field  at  the  outer  boundary,  it  becomes 
possible  to  reduce  the  number  of  unknowns  to  a  manageable  size  while  modeling  the 
physical  problem  as  correctly  as  possible. 

In  this  chapter,  a  general  background  on  the  finite  element  method  will  be  presented 
and  the  derivations  of  the  Bayliss,  Gunzburger,  and  Turkel  (BGT)  [17]  boundary  operators 
for  both  the  Helmholtz  and  Laplace  equations  will  be  reviewed. 

2.2  The  Finite  Element  Method 

In  the  finite  element  method  the  Helmholtz  or  Laplace  equation  is  solved  at  a  finite 
number  of  grid  points.  These  equations  are  discretized  through  the  use  of  a  weak 
variational  formulation.  In  this  section  only  Laplace's  equation  will  be  treated;  the 
extension  to  the  Helmholtz  equation  should  be  straightforward. 

Consider  the  N-conductor  configuration  shown  in  Figure  2.1.  The  potential,  u. 


satisfies  the  Laplace  equation  [1] 
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V-(eVu)  =  0  (2.1) 

where  £  =  e(x,y)  is  the  permittivity  of  the  dielectric.  The  following  boundary  conditions 
should  be  satisfied 

u  =  <5>j  on  r\ 
u  =  <p2  on  r2 

(2.2) 


u  =  <t>non  r„ 

Multiplying  (2.1)  by  a  testing  function  v  and  integrating  over  the  domain  of  the  problem  Q, 
we  obtain 

f  vV  (eVu)  ds  =  0  (2.3) 

Jn 


From  Green's  identity  we  have 


f  vV-(eVu)  ds  =  —  f  eVu'Vv  ds  +  f  ve  ^  dt 
Jn  Jsi  J  r0 


(2.4) 


Inserting  the  above  in  (2.3),  we  obtain 


f  eVu-Vvds-  f  ve^dt=0 
Jn  Jr  on 


(2.5) 
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Only  the  first  term  of  Equation  (2.5)  is  treated  in  this  section.  The  second  term  will  be 
examined  when  the  absorbing  boundary  condition  is  discussed. 

To  construct  an  approximate  solution,  the  domain  Cl  is  broken  into  a  number  of 
triangles.  Over  each  triangle  the  unknown  potential  is  expressed  in  terms  of  a  set  of 
approximating  basis  functions 

U(5i  »  S2  •  §3>  =  X  Uijk  *  S2  .  §3>  (2.6) 

ijk 

The  arguments  £1,  £2>  and  ^3  are  called  simplex  or  area  coordinates  and  can  be  related  to 
the  Cartesian  coordinates  by 


x  =  Si  xl+S2x2  +  S3  x3 

(2.7) 

y  =  Si  yi  +  S2  y2  +  S3  y3 

(2.8) 

where  xk  and  yk  are  the  coordinates  of  the  k*  vertex  shown  in  Figure  2.2. 

Note  that  each  approximating  function  is  equal  to  one  at  its  corresponding  node  and 
zero  at  all  other  nodes.  It  suffices  to  impose  continuity  of  the  functions  at  triangle  vertices 
to  guarantee  the  continuity  of  the  potential  across  all  of  the  triangle  edges.  It  should  also  be 
noted  that  the  simplex  coordinates  are  purely  local  in  nature  which  makes  the  analysis 
independent  of  the  position  of  the  triangular  elements  in  the  Cartesian  coordinate  system. 

In  order  to  determine  the  basis  function  otjjk,  we  need  to  first  define  a  family  of 
auxiliary  polynomials  R  of  degree  n.  Keeping  in  mind  that  ctijk  must  take  the  value  of 
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unity  at  node  ijk,  and  zero  at  all  other  nodes,  a  proper  choice  of  the  auxiliary  polynomials 
will  be 


k=0 


k=0 


R0(n^)  =1 


(2.9) 


The  polynomial  has  exactly  m  equispaced  zeros  at  ^  =  0,  1/n, (m-l)/n,  all  of  which  lie 
to  the  left  of  §  =  m/n,  and  takes  on  the  value  of  unity  at  1=  m/n.  A  family  of  interpolation 
functions  that  are  suitable  for  two-dimensional  problems  is  defined  by 


(Xjjk  =  RjCn.^jRjGi.l^RkfoJh)  i  +  j  +  k  =n  (2. 10) 

Clearly,  the  constructed  interpolation  functions  are  independent  of  the  global  coordinate 
system  because  they  are  defined  in  terms  of  simplex  coordinates.  Therefore,  the  creation  of 
the  finite  element  matrices  can  be  done  using  the  simplex  coordinates.  The  filling  of  the 
finite  element  matrices  can  then  be  placed  in  subroutines  that  can  be  called  for  each 
triangular  element  regardless  of  its  position  in  the  global  coordinate  system. 

Expressing  the  unknown  potential  u  as  in  Equation  (2.6),  the  first  term  of  equation 
(2.5)  will  take  the  form 

3-«X*S,  (2.11) 


where 


>ij  =  J  Va-Vaj 


ds 
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(2.12) 


The  element  matrix  S  may  be  written  in  detail  as 


K9cti  9a:  3a:  3a; 

irir  +  iir^Jds 


(2.13) 


which  may  be  expressed  in  terms  of  only  simplex  coordinates  and  Cartesian  coordinates  of 
the  vertices  of  each  triangle 


4A 


n  m 


=Sc°‘e*j  ( 


cncm)J* 

t 

3a{  3ctj 

ds 

da{ 

Y  3ctj 

9ctj  >  ds 

^k+i 

^k-i . 

3^k-i  J  2A 

(2.14) 


where  the  subscripts  n  and  m  progress  modulo  3,  A  is  the  area  of  the  triangle,  and 

bn=yn+i-yn-i  (2-15) 

c„  =  xn_i-xn+1  (2.16) 

0k  denotes  the  included  angle  at  vertex  k  and  it  is  given  by 


cot0k  =  - ^  (bnbm  +  cncm)  n*m 


(2.17) 


The  right-hand  side  of  Equation  (2.14)  is  dimensionless  and  depends  only  on  the  simplex 
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coordinates.  The  matrix  Sij  can  be  evaluated  for  a  triangle  of  any  size  or  shape  which  can 
then  be  stored  for  permanent  reference.  All  integrations  and  differentiations  can  be  carried 
out  in  a  universal  manner  valid  for  any  triangle  due  to  the  use  of  simplex  coordinates  [1], 

Each  element  matrix  is  then  calculated  and  added  to  the  global  system  of  equations 
on  an  element  by  element  basis  using  the  connectivity  matrix  which  has  information  on  the 
global  and  the  local  node  numberings.  If  two  nodes  numbered  1  and  i  are  not  in  adjaceru 
elements,  the  entry  (1,  i)  in  the  global  matrix  would  be  zero.  This  elegant  feature  of  the 
finite  element  method  leads  to  a  very  sparse  system  of  equations  that  can  be  easily  handled 
using  special  techniques.  Typically,  a  finite  element  problem  would  have  a  large  number  of 
unknowns  but  since  the  resulting  system  of  equations  is  very  sparse  it  can  be  solved 
efficiently.  The  truncation  from  an  infinite  region  to  a  finite  region  may  cause  a  significant 
error  in  the  solution.  The  best  known  method  of  dealing  with  this  problem  is  to  use  an 
absorbing  boundary  condition  at  the  outer  boundary.  In  the  following  section  a  review  of 
the  absorbing  boundary  conditions  will  be  given. 

2.3  The  Absorbing  and  Asymptotic  Boundary  Conditions 

When  using  FEM  to  solve  open  region  problems,  one  needs  to  introduce  an 
artificial  outer  boundary  in  order  to  limit  the  number  of  unknowns.  Since  the  outside 
region  is  unbounded,  one  needs  to  impose  the  corresponding  behavior  at  infinity  on  the 
finite  distance  boundary  and  obtain  an  accurate  solution  in  the  interior  region  [29].  Such 
boundary  conditions  that  dictate  the  behavior  at  infinity  are  commonly  referred  to  as 
asymptotic  or  absorbing  boundary  conditions  (ABC).  In  what  follows,  the  derivations  of 
the  Bayliss,  Gunzburger,  and  Turkel  (BGT)  [17]  boundary  operators  for  both  the 
Helmholtz  and  Laplace  equations  are  reviewed. 
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2.3.1  The  BGT  operator  for  the  Helmholtz  equation 

Consider  the  region  £2,  bounded  with  the  contour  T \,  shown  in  Figure  2.3.  Let  the 
exterior  region  to  Q  be  Hx  [29].  This  problem  is  equivalent  to 


V“u  +  k2u  =0  in  ftp 

(2.18) 

9u  r 

(2.19) 

u  satisfies  a  radiation  condition 


(2.20) 


where  u  is  the  scattered  field  and  g  is  the  contribution  from  the  incident  field.  In  the  far 
region  the  solution  of  (2.18)  takes  an  asymptotic  form 


u  = 


e-jkp 

Vp 


ao(<t>)- 


a^)  a2(<|>) 


eJ'kp 

Vp 


b0(<|))+ 


bj(4>)  t^Gf)) 


(2.21) 


The  first  term  of  Equation  (2.21)  designates  the  outgoing  waves  and  the  second  part 
designates  the  incoming  waves.  In  the  far  region  there  are  only  outgoing  waves;  therefore, 
the  second  term  of  Equation  (2.21)  is  not  physically  meaningful  and  only  the  first  term 
should  be  kept.  In  other  words,  the  scattered  field  in  the  far  region  should  behave  as 


u  = 


e*Y  ... 


(2.22) 


If  Up  designates  the  radial  derivative  of  u,  we  have  from  the  above  equation. 
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Incident  wave 


Figure  2.3.  Two-dimensional  scattering  problem. 
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e'jkp  ( 

Up+jku  =  — — 3^-^ao(4>)-»- 


a^tjO 
P  + 


a2(4>)  ^ 


c-tof&iW  2a 2«»)  ^ 

VF  l  p2  +  p3  +  ‘"J 


(2.23) 


Clearly,  Equation  (2.23)  is  equivalent  to 


(2.24) 


As  p  goes  to  infinity,  the  right-hand  side  of  Equation  (2.24)  goes  to  zero,  which  is 
equivalent  to  imposing  the  Sommerfeld  radiation  condition  on  the  field  u.  Examining  the 
first  term  of  Equation  (2.23),  we  note  that  it  is  equal  to  -u/2p.  Consequently,  we  can 
obtain  a  higher-order  boundary  condition  by  writing 

up+jta  +  ^-  =  °(^]  (2'25 


This  type  of  analysis  was  first  carried  out  by  Bayliss  and  Turkel  [21]  as  well  as  Engquist 
and  Majda  [16]  in  connection  with  time -dependent  problems.  For  time  harmonic  problems, 
their  results  reduce  to  Equation  (2.25).  Bayliss,  Gunzburger,  and  Turkel  [17]  have 
generalized  these  results  by  writing 


Bi*^+jk+2? 


(2.26) 


or,  equivalently, 


,,-jkp 


a  j  (<]>)  2a2($) 

v  +~~3 
v  p  p 


j 


(2.27) 
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They  were  able  to  obtain  a  higher-order  operator  by  letting  v  =  Biu  and  observing  that 


3v  5  J  1 

3p+*'3pv-°b* 


U1 

B^u = ( J? + jk + + ik + ^)u = 


(2.28) 


Funhermore,  they  were  able  to  generalize  the  procedure  by  deriving  an  m1*1  operator ,  Bm, 


such  that 


BmU  p2m+l/2 


(2.29) 


fa  H)  ,k 

■5—  + - +  jk 

V3p  p  J 


(2.30) 


The  original  problem  has  now  been  slightly  modified.  That  is,  the  infinite  exterior  region, 
Qx,  is  now  truncated  and  bounded  by  a  circular  contour  T2  where  the  m*  operator  can  be 
applied.  Obviously,  one  can  now  use  a  PDE  technique  to  solve  the  approximate  problem 
provided  the  region  bounded  by  T2  is  not  so  large  as  to  require  an  unmanageable  number  of 
unknowns.  The  approximate  problem,  shown  in  Figure  2.4,  is  equivalent  to 


V2u  +  k2u  =  0  in  Qt 


(2.31) 
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9u 

^  =  gon 


(2.32) 


Bmu  =  0  on  r2  (2.33) 

Note  that  the  m*  order  operator  has  an  mth  order  radial  derivative.  By  invoking  the 
Helmholtz  equation  we  can  trade  the  mlh  order  radial  derivatives  for  the  first-order  radial 
derivatives  plus  derivatives  in  the  tangential  directions.  Note  also  that  there  is  an  error 
introduced  by  the  approximate  problem  of  the  order  0(1/R2m+1/2),  where  R  is  the  distance 
from  the  outer  boundary  to  the  origin.  In  two  dimensions,  it  is  impossible  to  know  the 
corresponding  error  in  ilu-uaproH.  where  uapro  is  the  solution  of  the  approximate  problem. 
In  three  dimensions,  however,  Bayliss,  Gunzburger,  and  Turkel  were  able  to  prove,  for 
m=l  and  m=2,  the  following  theorem 


(2.34) 


where  Udif  =  u-uapr0,  C  depends  on  k  and  the  outer  boundary,  and  the  surface  norm  is 
defined  as 


ds 


(2.35) 


Hence,  the  error  in  the  solution  caused  by  the  truncation  of  the  infinite  domain  is  the  same 
on  the  artificial  boundary  as  on  the  interior  boundary  and  is  inversely  proportional  to 
l/rm+1. 
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2.3.2  The  BGT  operator  for  the  Laplace  equation 

Again  consider  the  problem  described  in  Figure  2.3.  This  time  we  would  like  to 
consider  Laplace's  equation  which  means  that  there  is  neither  an  incident  nor  a  scattei  d 
wave.  The  problem  at  hand  is  equivalent  to 

V2u  =  OinOr  (2.36) 

usgonTj  (2.37) 

u  =  logp  as  p— ><»  (2.38) 

where  u  is  the  electrostatic  potential. 

The  general  solution  of  this  problem  can  written  as 

u(p,4>)  =  logp  +  ^  ~  cosn<{>  (2.39) 

n=0  P 


or  more  explicitly  as 


(2.40) 


We  will  now  derive  a  set  of  boundary  condition  operators  that  can  be  applied  ~>n  the 
artificial  outer  boundary.  If  we  let  v  =  u  -  logp  -  aq,  we  can  see  that 
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9v 

= 


3j  2a2 '  3a3 

- -cos<}> - — cos2<{t - —  cos3<j>  -  ... 

P  P  P 


(2.41) 


It  can  be  seen  that 


3v  v  i  (  2a5 

3p  +  —  =  — 5\a2cos2<!>  +  —  cos30  +  ...J  (2.42) 


Thus, 


We  then  define  the  first-order  operator  to  be 

_  a  i 
B‘$  ¥+p 


(2.43) 


(2.44) 


The  second-order  operator  can  be  obtained  by  letting  w  =  Biv  and  observing  that 


9w 

ap 


+ 


3w 

P 


2a3 

=  — -cos3<j)  +  ... 

P 


(2.45) 


It  can  readily 


B2  = 


be  verified  that 


(2.46) 


The  process  can  be  repeated  to  obtain  the  mth  order  operator  which  can  be  written  as 
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(2.47) 


The  infinite  exterior  region  Qj  is  now  truncated  and  bounded  by  a  circular  contour  F2 
where  the  111th  operator  can  be  applied.  The  approximate  problem  described  in  the  truncated 
region,  shown  in  Figure  2.4,  is  equivalent  to 


V2u  =  0  in  QT 

(2.48) 

u  =  g  on 

(2.49) 

Bmu  =0  on  r2 

(2.50) 

In  the  region  Qj  one  can  now  use  a  PDE  technique  to  solve  the  approximate  problem, 
provided  that  the  region  bounded  by  T2  is  not  so  large  as  to  require  an  unmanageable 
number  of  unknowns. 


2.4  Conclusions 


In  this  chapter,  a  brief  overview  of  the  finite  element  method  and  its  implementation 
were  presented.  Starting  from  the  asymptotic  representation  of  the  solution  for  the 
Helmholtz  and  Laplace  equations,  the  Bayliss,  Gunzburger,  and  Turkel  (BGT)  mlh  order 
absorbing  and  asymptotic  boundary  condition  operators  were  derived  for  both 
electromagnetic  scattering  and  electrostatic  problems.  In  Chapter  3,  the  BGT  operator  for 
the  Helmholtz  equation  will  be  used  to  derive  an  absorbing  boundary  condition  operator 
that  can  be  applied  on  an  arbitrary  boundary.  The  resulting  operator  will  then  be  used  to 
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study  the  scattering  from  some  long  scatterers  where  the  outer  boundary  will  be  arbitrary. 
In  Chapter  4,  a  similar  operator  for  electrostatic  problems,  based  on  the  BGT  operator 
associated  with  Laplace's  equation,  will  be  derived. 
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CHAPTER  3 

AN  ABSORBING  BOUNDARY  CONDITION 
FOR  AN  ARBITRARY  OUTER  BOUNDARY 

3.1  Introduction 

The  finite  element  method  is  very  appealing  for  solving  open  region  scattering 
problems  because  of  its  simplicity  in  modeling  complex-shaped  structures  and 
inhomogeneous  dielectrics.  The  local  or  the  absorbing  boundary  condition  makes  the  FEM 
even  more  powerful  because  it  preserves  the  sparsity  of  the  matrix.  However,  in  the  past, 
there  has  been  a  major  limitation  associated  with  the  ABC  as  applied  to  scattering  problems 
[32],  That  is,  most  of  the  attempts  to  solve  the  open  region  scattering  problem  using  the 
local-type  boundary  operators  have  been  limited  to  separable  types  of  geometry  for  the 
outer  boundary.  This  may  be  attributable  to  the  form  in  which  these  operators  were  cast  in 
their  original  versions.  However,  for  a  class  of  scatterers  that  are  long  and  slender,  e.g.,  a 
strip  or  an  airfoil,  truncating  the  open  region  with  a  circular  outer  boundary  requires 
solving  for  the  field  over  a  very  large  mesh  region  surrounding  the  scatterer,  which,  in 
turn,  requires  the  solution  of  a  large  matrix.  In  such  situations,  the  FEM  becomes  less 
efficient  than  the  method  of  moments.  However,  if  the  outer  boundary  is  made  to  conform 
to  the  geometry  of  the  scatterer,  a  significant  improvement  in  computation  time  and  storage 
can  be  achieved. 

In  this  chapter,  we  show  how  we  can  use  an  outer  boundary  of  an  arbitrary  shape 
that  can  be  made  to  conform  to  the  geometry  of  the  scatterer  in  order  to  minimize  the  size  of 
the  solution  region.  As  will  be  shown  below,  this  requires  transformation  of  the  boundary 
operator  into  a  form  that  uses  a  local  coordinate  system.  Once  this  is  done,  the  operator  can 
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be  applied  on  a  per  element  basis  in  the  FEM  formulation.  To  illustrate  the  application  of 
the  newly  transformed  boundary  operator,  we  present  the  results  for  several  p.e.c. 
scatterers. 

3.2  Derivation  of  the  Boundary  Condition 

In  the  previous  chapter,  a  brief  discussion  of  the  finite  element  method  and  the 
absorbing  boundary  condition  was  presented.  A  detailed  explanation  of  the  implementation 
of  the  absorbing  boundary  condition  in  the  finite  element  formulation  for  a  two-dimensional 
scatterer  enclosed  by  a  circular  outer  boundary  was  given  in  [24].  For  a  complete 
coverage  of  the  finite  element  method  and  absorbing  boundary  condition,  the  reader  is 
referred  to  the  work  of  various  authors  that  appeared  in  [8]-[32].  In  this  section,  we 
develop  a  finite  element  scheme  for  solving  the  problem  of  scattering  by  objects  of  arbitrary 
shape. 

3.2.1  Formulation 

Consider  an  arbitrarily-shaped  perfectly  conducting  scatterer  whose  exterior  region, 
Q,  is  bounded  by  the  contour  Ti,  as  shown  in  Figure  3.1.  For  a  TM  or  TE  polarized 
incident  wave,  the  scattered  field,  u,  satisfies  the  wave  equation 

(V2  +  k2)u  =  0  (3.1) 

To  obtain  the  variational  expression  for  this  equation,  we  multiply  (3.1)  by  a  testing 
function,  v,  and  integrate  over  the  domain  of  the  problem,  Q,  to  obtain 
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Incident  wave 


Figure  3.1.  Geometry  for  the  finite-mathematics  approach  to  the  scattering 
problem  where  the  outer  boundary  is  arbitrary. 
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u  +  k2vu)  ds  =  0 


(3.2) 


Using  the  Green’s  identity  to  transfer  the  differentiation  from  the  unknown  function  u  to 
the  testing  function  v,  we  obtain 

f  vV2uds  =  -f  Vu-Vvds+f  v^-dl+f  v^-dl  (3.3) 

Jn  Jrt  °n  Jr2  °n 

In  the  case  of  TM  polarization,  where  the  total  field  at  the  p.e.c.  scatterer’s  surface  is  zero, 
there  will  not  be  any  contribution  from  the  boundary  integral  on  Ti,  because  the  field 
values  are  specified  on  the  surface  of  the  scatterer. 

Substituting  (3.3)  in  (3.2) ,  we  have 

J  (Vu-Vv  -  k2vu)  ds  =  f  v-|~  dl  +  J  v^-  dl  (3.4) 

*i  2 

The  second  term  of  the  right-hand  side  of  Equation  (3.4)  involves  an  integral  over  the  outer 
boundary,  1*2,  and  the  normal  derivative  of  u  appears  in  this  integrand.  Obviously,  the 
absorbing  boundary  condition  has  to  be  applied  on  the  outer  boundary  T2-  Since  the 
normal  derivative  of  u  appears  in  the  boundary  integral  contribution  ,  the  absorbing 
boundary  condition  has  to  be  imposed  on  that  quantity.  Hence,  for  our  purposes,  it  is 
more  desirable  to  find  an  asymptotic  representation  for  the  normal  derivative  of  u  rather 
than  make  direct  use  of  the  BGT  operator  B2  as  given  in  (2.28).  For  a  circu’ar  outer 
boundary,  the  normal  derivative  is  simply  the  radial  one.  Using  the  BGT  operator  B2  as 
given  in  (2.28)  in  conjunction  with  (3.1),  we  obtain  an  asymptotic  representation  for  the 
radial  derivative  that  reads 
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u-,  =  a(p)  u  +  (3(p)  u, 


00 


(3.5) 


where  u  is  the  scattered  field,  is  the  second-order  angular  derivative,  and  a(p)  and 
P(p)  are  given  by  the  following  asymptotic  series  : 


,  v  1  j  1 

a(p)  =  -jk  -  -r- - r  + 


2P  8kp2  8k2p3 


(3.6) 


p(p)  -  +  1 


2kp2  2k2p3 


(3.7) 


In  a  previous  work  [24],  T2  was  chosen  to  be  a  circle  enclosing  the  scatterer. 
Thus,  the  normal  derivative  of  u  was  simply  its  radial  derivative  up.  Our  next  task  is  to 
transform  the  absorbing  boundary  operator,  as  given  in  (3.5),  into  a  form  suitable  for 
arbitrary  boundaries.  In  view  of  (3.4),  it  is  natural  to  attempt  to  write  the  new  operator  in 
terms  of  the  normal  derivative  of  the  field  on  the  boundary.  Since  the  region  of  solution 
will  be  discretized  into  finite  elements,  the  transformed  operator  will  be  made  local,  i.e.,  the 
normal  derivative  on  the  outer  boundary  will  depend  locally  on  the  coordinates  of  each 
element.  This  will  be  demonstrated  below  when  we  derive  the  expression  for  the 
transformed  operator. 

Consider  the  triangle  shown  in  Figure  3.2.  The  edge  1-2  resides  on  the  outer 
contour  T2.  Let  n-t  represent  a  new  coordinate  system  with  the  origin  at  node  1  with 
coordinates  (xo,  yo).  The  relationship  between  the  cartesian  coordinate  system  and  the  new 
local  coordinate  system  can  be  obtained  via  the  following  equations  that  pertain  to  a 
translation  and  rotation  of  the  coordinate  system: 


n  =  (x-x0)  sin0o  -  (y-y0)cos90 


(3.8) 


t  =  (x-xo)cos0o  +  (y-yo)sin0o 

(3.9) 

x  =  nsin90  +  tcos90  +  Xq 

(3.10) 

y  =  -ncos90  +  tsin0o  +y0 

(3.11) 

where  9o  is  the  angle  the  tangent  vector  makes  with  the  positive  x-axis. 

Our  aim  is  to  obtain  an  expression  for  the  normal  derivative  un  on  the  t-axis  in 
terms  of  tangential  derivatives.  Using  the  chain  rule  on  un,  we  have 

3u  _  3u  9p  9u_90 
3n  ~  9p  3n  +  9$  9n 


where  p  and  p  are  the  cylindrical  coordinates  of  a  point  on  the  outer  boundary  I~2. 

To  this  end,  the  radial  derivative  up  is  given  by  (3.5).  Expressions  for  the 

remaining  partial  derivatives  are  then  obtained: 

3n  "  3x  3n  +  9y  3n  1 

8<t>  _  9<J)  3x  9<t>  3y 
3n  ~  9x  ^n  9y  3n 


For  the  triangular  element  edge  residing  on  the  t-axis,  the’normal  component  in  the  n-t 
coordinate  system  is  zero.  Consequently,  Equations  (3.10)  and  (3.1 1)  become 


32 


X  =  tcos0o  +  x0 

(3.15) 

y  =  tsin0o  +  y0 

(3.16) 

Thus,  we  obtain: 

9p  tcos  0O  +  x0 

9x  /  2  2  2 

y  t  +  xq  +  y0  +  2t  (xocos0o  +  yosin0o) 

(3.17) 

dp  tsin0o  +  y0 

ay  /  2  2  2 

t  +  Xo  +  yo  +  2t  (xqCOS0q  4-  yosin0o) 

(3.18) 

9<j>  tsin0o  +  y0 

t2  +  x^  +  yo  +  2t  (xocos0o  +  yosin0o) 

(3.19) 

tcos0o  +  x0 

t2  +  +  yo  +  2t  (xocos0o  +  yosin0o) 

(3.201 

In  order  to  obtain  an  expression  for  the  angular  derivative  u<j). 

we  use  the  following 

approximation: 

9u  9u  dt 

90  9t  90 

(3.21) 

Using  the  chain  rule,  we  obtain 

|^  =  *o  sin0o  -  y0  cos0o 

(3.22) 
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Consequently,  the  first-order  angular  and  tangential  derivatives  are  related  by 


Ij-  -  (x°  sin60  -  y0  cos90) 


(3.23) 


and  the  second-order  angular  and  tangential  derivatives  by 


2  2 

=  ~r  (xo  sin0o  -yo  cose0)2 

a2<h  a  t 


(3.24) 


Using  the  expressions  of  the  various  partial  derivatives  obtained  so  far  along  with 
Equations  (3.5)-(3.7),  we  arrive  at  the  final  form  of  the  normal  derivative  un 


|^  =  au+y  ut  +  P  uu 


(3.25) 


where 


a  =  (x0  Sin90  -  y0  cos60)f-“  ~  +  -j-j 

V  K  2p  8kp  8k  p 


(3.26) 


a  2  2 

_  -t  (x0  sin0o  -y0  cos0o)  +-jsin20o  (y0  -x0)  +xQy0  cos20o 


(3.27) 


—  3f  j  1 

p  =  (x0  sin60  -y0  cos0o) - -  +  — — 

\  2kp  2k  p 


(3.28) 


where  p  is  given  by 


V2  2  2 

t  +xo +  yo +  2t  (x0coseo  +yosin0o^ 


(3.29) 
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where 

C  =  (xosin0o-yocos0o)3(t  +  xocos0o+yosin0o/  — - - (3.32) 

l 2kp  k  p 

The  form  given  in  (31)  is  well-suited  for  numerical  implementation  for  any  arbitrary 
scatterer  enclosed  by  any  arbitrarily-shaped  outer  boundary  I>. 

3.2.2  Radar  cross-section  calculation 


The  radar  cross-section  is  an  important  quantity  that  needs  to  be  computed  for 
scattering  problems.  It  is  defined  in  [33]  as  the  area  for  which  the  incident  wave  contains 
sufficient  power  to  produce,  by  omnidirectional  radiation,  the  same  back-scattered  power 
density. 
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3.2.2. 1  TM  case 


For  TM  incidence,  the  two-dimensional  radar  cross-section  is 


CTTM^-0  =  Lip^2nP 


e;(P,o>)12 

Ez(0,0)|2 


(3.33) 


where  (p,<J>)  are  polar  coordinates.  For  a  plane  wave  incident  of  the  form 


The  scattered  field  is  equal  to 


9F  3F 

E’(x,y)  =  -jfcnA2--5^  +  1;r 


where 

Az(x,y)  =  f  Jj2(x’,  y')  Hg2)(kR)  dx'dy' 
F(x,y)  =  |J  K(x’,  y')  H^CkR)  dx'dy' 


and 


„  /  ,  2 ,  ,  A  ^  (x')2  +  (y')2 

R  =  p<y  /  1 — (x  cos<J)  +  y  sinp)  + - - - 

P 


(3.34) 


(3.35) 


(3.36) 

(3.37) 


(3.38) 


In  the  far  region,  the  third  term  under  the  radical  is  negligible.  Thus  using  the 
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approximation 


V  1+e  =  1  +  j 


R  takes  the  form 


R  =  p  -  x'costj)  -  y’sin4>  as  p— 


Using  the  asymptotic  form  of  the  Hankel  function 


r  T(2)  <•  , 
H0  (a) 


2j  -ja 

,=r-  e  as  a—>°° 

ria 


(3.39) 


(3.40) 


(3.41) 


in  Equations  (3.35M3.37),  we  obtain  an  expression  for  the  radar  cross-section 

OTM(0,  4>ij,c)  =  j  I J  J 0UZ  +  Kxsin<J>  -  Kycos<J>)ejk(x'cos0+y'sino)dx'dy'  I  2  (3.42) 


For  a  circular  outer  boundary,  an  FFT  algorithm  is  usually  used  to  calculate  the  radar  cross- 
section.  However,  for  an  arbitrary  outer  boundary  it  is  not  possible  to  use  the  FFT 
method,  and  (3.42)  needs  to  be  used  to  approximate  the  radar  cross-section.  Since  we  are 
dealing  with  TM  scattering,  we  are  basically  solving  for  the  z-component  of  the  scattered 
field.  In  order  to  find  the  electric  current  that  appears  in  (3.42),  it  would  seem  reasonable 
to  use  the  finite  difference  method  to  obtain  the  numerical  derivative  of  the  scattered  electric 
field.  However,  we  found  that  taking  a  numerical  derivative  is  very  unstable,  especially  in 
the  region  around  the  scatterer  where  the  near  field  has  many  oscillations.  To  circumvent 
this  problem,  we  have  used  the  absorbing  condition  operator  given  by  (3.25)  to  find  the 
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normal  derivative  of  the  scattered  field  on  the  outer  boundary  from  which  the  equivalent 
electric  current  is  calculated.  It  should  be  understood  that  the  use  of  the  normal  derivative 
expression  given  by  (3.25)  is  also  an  approximation,  but  this  method  is  found  to  be  more 
stable  than  taking  a  numerical  derivative  via  the  finite  difference  method. 

3.2.2.2  TE  case 

For  TE  incidence,  the  two-dimensional  radar  cross-section  is 


®TE^’^inP  =  Li  m  .  j 

1 E  mc  o— **> 


H*(p,p)  | 


P^~  H;(0,0) 


(3.43) 


For  a  plane  wave  incident  of  the  form 


H*(x,y)  =  e^k<xcos^"‘+  ysinC".) 


(3.44) 


the  scattered  magnetic  field  is  given  by 


_  9Ay  9Ax  .  k 


H:(x,y)  =  -5-i  -  -y-i  -  j-  F7 
1  }  ox  dy  J  T|  z 


(3.45) 


Using  the  same  far-field  approximation  as  the  one  used  earlier  for  the  TM  case,  we  obtain 


dx  dy  4  v  FLkp  J  J 


f  f  (J  cos$-J,sin({0e^k<*cos4+y  sind\ix'dy'  (3.46) 


(3.47) 
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Combining  (3.43)-(3.47)  gives  the  radar  cross-section 


-jiJK 


Kz 

J  sin<t>  -  Jcos<S>  H - 


ejk(x-cosO+ysinO)dx,dyl 


(3.48) 


In  the  case  of  a  p.e.c.  scatterer,  the  equivalent  sources  on  the  surface  of  the  scatterer  are 
purely  electric  and  hence  the  radar  cross-section  simplifies  to 

=  J I J  J  (Jxsin<|)-Jycos(j))ejk(x'cos0+y'sinOdx*dy’  1 2  (3.49) 

For  TE  scattering,  one  solves  for  the  z-component  of  the  scattered  magnetic  field  which  is 
directly  proportional  to  the  electric  current  on  the  scatterer.  Therefore,  the  implementation 
of  (3.49)  to  calculate  the  radar  cross-section  for  a  p.e.c.  scatterer  does  not  pose  any 
numerical  difficulty. 

3.3  Numerical  Results 

As  we  indicated  earlier,  the  form  given  in  (3.31)  is  well-suited  for  numerical 
implementation  for  any  arbitrary  scatterer  enclosed  by  any  arbitrarily-shaped  outer 
boundary  r2.  Equation  (3.31)  was  used  to  investigate  the  problem  of  scattering  by  several 
p.  e.  c.  scatterers  using  some  elongated  and  conformable  outer  boundaries.  For  mesh 
generation,  Patran  was  used  for  all  of  the  cases  considered.  Although  (3.31)  is  valid  for 
bn y  order  triangular  elements,  only  first-order  elements  were  used.  For  the  TM  case,  a 
mesh  density  of  13  to  15  nodes  per  wavelength  was  used.  Whereas  for  the  TE  case,  to 
achieve  a  reasonable  accuracy,  a  mesh  density  of  15  to  17  nodes  per  wavelength  was 
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required.  It  was  also  noticed  that  denser  meshes  are  required  for  larger  scatterers. 

3.3.1  TM  incidence  (E«  =  e*Jkx) 

The  finite  element  solution  gives  the  z-component  of  the  scattered  electric  field.  As 
indicated  earlier,  there  is  an  added  error  in  the  radar  cross-section  computation  due  to  the 
difficulties  encountered  in  approximating  the  normal  derivative  at  the  outer  boundary  and, 
hence,  in  implementing  (3.42).  The  first  scatterer  considered  was  a  4X  strip,  shown  in 
Figure  3.3,  enclosed  by  an  elongated  boundary  IX  away  from  the  end  of  the  strip  in  the  x- 
direction  and  IX  away  in  the  y-direction.  The  second  scatterer  was  a  2X  by  IX  wedge,  for 
which  the  outer  boundary  was  chosen  to  be  as  conformable  as  possible  to  the  surface  of  the 
scatterer,  shown  in  Figure  3.4.  The  radar  cross-section  was  calculated  using  the  present 
method  and  the  method  of  moments,  and  the  results  are  shown  in  Figures  3.5  and  3.6.  In 
both  cases,  the  solutions  obtained  using  the  present  method  agree  favorably  with  those 
derived  via  the  method  of  moments.  The  third  scatterer  considered  was  a  9X  strip  enclosed 
by  an  elongated  boundary  IX  away  from  the  end  of  the  strip  in  the  x-direction  and  IX  away 
in  the  y-direction,  shown  in  Figure  3.7.  The  near  field  was  calculated  on  the  outer 
boundary  using  the  present  method  and  the  method  of  moments  and  plotted  versus  the 
angle  <J)  where  <J>  is  as  shown  in  Figures  3.2  and  3.7  (Figure  3.8).  When  the  outer 
boundary  was  extended  to  2X  away  in  the  y-  direction,  as  shown  in  Figure  3.7b,  the 
solution  obtained  using  the  present  method  agreed  favorably  with  that  computed  with  the 
method  of  moments  (Figure  3.9).  The  fourth  scatterer  considered  was  a  8X  by  4X  p.e.c. 
wedge,  and  the  outer  boundary  was  chosen  to  be  conformable  to  the  surface  of  the  scatterer 
as  shown  in  Figure  3.10.  As  Figure  3.1 1  indicates,  the  near  field  on  the  outer  boundary 
given  by  the  finite  element  method  agrees  reasonably  well  with  the  method  of  moments' 
solution. 


Figure  3.3.  A  4X  strip  enclosed  with  an  elongated 
the  x-  and  y-directions. 


boundary  I~2  IX.  awa 
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Figure  3.4.  A  2X  by  IX  wedge  enclosed  by  an  outer  boundary  T2 
having  a  general  shape  similar  to  that  of  the  scatterer. 


0  20  40  60  80  100  120  140 

Angle  in  degrees 


Figure  3.5.  Radar  cross  section  of  a  4X  strip  illuminated  by 
a  TM  incident  wave  and  enclosed  by 

conformable  outer  boundary  1  X  away  from  the 
surface  of  the  scatterer  as  shown  in  Figure  3.3. 
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Figure  3.6.  Radar  cross  section  of  2\  by  IX  wedge 
illuminated  by  a  TM  incident  wave  and 
enclosed  by  a  conformable  outer  boundary 

0.3X  away  from  the  surface  of  the  scatterer  as 
shown  in  Figure  3.4. 


Figure  3.7.  A  9X  strip  enclosed  with  an  elongated  outer  boundary  T2. 

(a)  T2  is  IX  away  from  the  surface  of  the  scatterer  in  the  x-  and 
y-directions. 

(b;  T2  is  IX  away  from  the  surface  of  the  scatterer  in  the  x-direction 
and  2X  in  the  y-direction. 
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(a) 


0  18  36  54  72  90  108  126  144  162  180 


Angle  in  degrees 
(b) 

Figure  3.8.  Near  field  on  the  outer  boundary  of  a  9  A.  strip  illuminated  by  a 
TM  incident  wave  and  enclosed  by  an  elongated  boundary  as 
shown  in  Figure  3.7(a). 

a)  Real  part  of  near  field. 

b)  Imaginary  part  of  near  Field. 


46 


Angle  in  degrees 
(a) 


Angle  in  degrees 
(b) 

Figure  3.9.  Near  field  on  the  outer  boundary  of  a  9  X  strip  illuminated  by  a 
TM  incident  wave  and  enclosed  by  an  elongated  boundary  as 
shown  in  Figure  3.7(b). 

a)  Real  part  of  near  field. 

b)  Imaginary  part  of  near  field. 


Real  part  of  near-field 
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0  18  36  54  72  90  108  126  144  162  180 


Angle  in  degrees 

(a) 


Angle  in  degrees 
(b) 

Figure  3.11.  Near  field  on  the  outer  boundary  of  an  8X.  by  4X,  wedge 
illuminated  by  a  TM  incident  wave  and  enclosed  by  an 
elongated  boundary  as  shown  in  Figure  3.10. 

a)  Real  part  of  near  field. 

b)  Imaginary  part  of  near  field. 
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3.3.2  TE  incidence  (H»  =  e’jky) 

The  finite  element  solution  gives  the  z-component  of  the  scattered  magnetic  field. 
Since  the  electric  current  on  the  p.e.c.  scatterer  is  directly  proportional  to  the  magnetic  field, 
the  implementation  of  (3.49)  does  not  represent  any  numerical  difficulty.  The  wedge 
shown  in  Figure  3.4  was  illuminated  by  a  TE  wave.  For  this  case,  the  radar  cross-section 
has  been  computed.  As  is  evident  from  Figure  3.12,  our  results  compare  favorably  with 
those  obtained  via  the  method  of  moments. 

In  order  to  achieve  a  reasonable  accuracy  for  a  slender  scatterer  like  a  strip,  the 
outer  boundary  has  to  be  placed  a  distance  between  L/10  and  2L/10  in  the  x-direction  and 
between  2L/10  and  3L/10  in  the  y-direction  away  from  the  surface  of  the  scatterer  where  L 
is  the  total  leng'h  of  the  strip.  For  instance,  to  achieve  a  reasonable  accuracy  for  the  9X 
strip  shown  in  Figure  3.7b,  it  was  required  to  place  the  elongated  outer  boundary  13.  in  the 
x-diKCtion  and  23.  in  the  y-direction  away  from  the  surface  of  the  scatterer.  Using  a 
density  of  225  nodes/3.2  and  optimizing  the  mesh,  the  93.  strip  shown  in  Figure  3.7b 
required  a  mesh  size  of  8200  .nodes.  If  one  were  to  enclose  the  93.  strip  with  a  circular 
outer  boundary  having  a  radius  of  5.53.  and  using  the  same  mesh  density  as  the  one  used 
for  the  elongated  outer  boundary,  one  would  end  up  with  a  mesh  size  of  19100  nodes.  For 
a  sparse  matrix  solver  the  factorization  time  is  proportional  to  N2  while  the  matrix  fill  time 
is  proportional  to  N,  where  N  is  the  total  number  of  nodes.  Therefore,  the  use  of  an 
elongated  instead  of  a  circular  outer  boundary  resulted  in  a  reduction  by  a  factor  of  5.43  in 
the  factorization  time  and  a  factor  of  2.33  in  the  matrix  fill  time.  Furthermore,  the  storage 
requirement  is  also  reduced  by  at  least  a  factor  of  2.33.  Clearly,  there  is  a  distinct 
advantage  in  using  an  elongated  versus  a  circular  outer  boundary. 

As  the  numerical  results  indicate,  while  the  finite  element  yielded  acceptable  results 
for  all  the  strips  and  wedges  considered,  the  results  for  the  wedges  agreed  more  favorably 
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than  those  of  the  strips  with  the  method  of  moments  results.  This  can  be  explained  by  the 
fact  that  the  scattered  waves  are  purely  outgoing  only  in  the  region  outside  the  smallest 
circle  that  entirely  encloses  the  scatterer.  For  a  wedge,  where  the  outer  boundary  resembles 
more  closely  a  circular  one,  there  are  more  points  satisfying  the  above  criterion  (Figure 
3.13)  than  there  are  for  a  strip  where  the  outer  boundary  is  elongated  (Figure  3.14).  In 
fact,  the  presence  of  incoming  waves  is  more  pronounced  for  the  region  of  the  elongated 
boundary  enclosing  the  strip  which  is  close  to  the  origin. 
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Angle  in  degrees 

Figure  3.12.  Radar  cross  section  of  a  2k  by  IX  wedge  illuminated  by  a  TE 
incident  wave  and  enclosed  by  a  conformable  outer  boundary 
0.3Xaway  from  the  surface  of  the  scatterer  as  shown  in  Figure 
3.4. 


Smallest  circle 
enclosing  the  strip 


Figure  3.14.  Region  of  purely  outgoing  waves  for  a  strip. 
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3.4  Conclusions 

In  this  chapter  we  have  introduced  a  technique  for  mesh  truncation  using  an 
arbitrary  outer  boundary  and  have  derived  the  appropriate  absorbing  boundary  condition 
(ABC)  for  such  a  boundary.  An  important  feature  of  the  new  boundary  condition  is  that  it 
enables  one  to  use  a  boundary  that  conforms  to  the  geometry  of  the  scatterer.  The 
formulation  has  been  verified  by  comparing  the  results  with  those  obtained  for  the  method 
of  moments  for  scattering  from  a  4X  strip,  a  2X  by  IX  wedge,  a  9X  strip,  and  an  8X  by  4X 
wedge.  We  believe  that  this  boundary  condition  makes  it  practical  to  solve  the  arbitrarely- 
shaped  large-body  scattering  problems  using  FEM  by  reducing  the  number  of  mesh  points 
to  a  more  manageable  size  than  would  be  possible  with  a  circular  outer  boundary.  In 
Chapter  6,  an  absorbing  boundary  condition  that  takes  into  account  both  the  lower-  and 
higher-order  harmonics,  the  higher-order  absorbing  boundary  condition,  will  be  addressed. 
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CHAPTER  4 

AN  ASYMPTOTIC  BOUNDARY  CONDITION  FOR  QUASI-TEM 
ANALYSIS  OF  TWO-DIMENSIONAL  TRANSMISSION 
LINE  STRUCTURES 


4.1  Introduction 

Microwave  transmission  lines  have  been  investigated  by  many  researchers  who 
have  employed  a  variety  of  methods  to  study  the  problem  of  computing  the  characteristic 
impedance  and  propagation  constant  along  these  lines.  Some  of  these  techniques  include 
the  Fourier  transform  method  [34]-[35],  variational  method  [36]-[37],  spectral  domain 
method  [38]-[39],  Green's  function  technique  [40]-[46],  [51],  and  [54]  conformal 
mapping  [47]-[48],  boundary  element  method  [49]-[50],  and  Finite-element  method  [2], 
[3],  and  [7].  All  but  the  last  three  approaches  mentioned  above  are  limited  to  thin  strips 
and/or  to  structures  containing  dielectrics  with  planar  interfaces.  Although  the  finite 
element  method  method  (FEM)  is  very  general,  and  can  handle  any  arbitrary  configuration 
of  conductors  and  dielectrics,  it  must  deal  with  the  practical  problems  of  mesh  truncation 
and  the  need  for  a  large  number  of  mesh  nodes  when  applied  to  an  open  region  problem. 
One  approach  to  circumventing  this  difficulty  is  to  truncate  the  mesh  by  introducing  a 
fictitious  conducting  enclosure  [2]  and  [7].  This  approach  yields  satisfactory  results  only  if 
the  actual  field  decays  sufficiently  well  as  it  reaches  the  outer  boundary.  Typically,  this 
requires  one  to  move  the  outer  boundary  far  away  from  the  structure  in  order  to  achieve 
acceptable  accuracy  and,  this,  in  turn,  results  in  a  large  mesh.  An  alternative  approach  is  to 
use  "infinite”  elements  [3],  which  extend  to  infinity,  and  cover  the  region  outside  of  a 
fictitious  boundary  surrounding  the  structure.  Although  superior  to  the  artificial  p.e.c 
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boundary  method,  this  approach  nonetheless  has  its  own  drawbacks.  First,  the  infinite 
elements  require  special  care  during  the  filling  of  the  FEM  matrix.  Second,  one  needs  to 
assume  a  certain  asymptotic  behavior  of  the  field  within  the  infinite  elements,  and  this 
behavior  may  not  be  convenient  to  obtain. 

In  this  chapter,  we  introduce,  once  again,  the  concept  of  an  asymptotic  boundary 
condition  which  provides  us  with  an  efficient  means  for  dealing  with  the  open  region 
problems  in  the  quasi-static  regime.  This  asymptotic  boundary  condition  does  not  suffer 
from  the  complications  associated  with  the  infinite  elements,  and  yet  enables  us  to  bring  the 
outer  boundary  much  closer  to  the  structure  than  would  be  possible  with  the  p.e.c.  artificial 
boundary.  Furthermore,  unlike  many  of  the  available  ABCs  that  are  restricted  to  separable 
outer  boundaries,  the  one  presented  in  this  chapter  is  useful  for  an  arbitrarily-shaped  outer 
boundary.  We  will  demonstrate  the  versatility  of  this  new  asymptotic  boundary  condition 
by  considering  the  examples  of  one,  two,  and  six  conductor  microstrip  lines. 

4.2  Derivation  of  the  Asymptotic  Boundary  Condition 

Figure  4.1  depicts  the  geometry  of  an  open  region  problem  consisting  of  N 
arbitrarily-shaped  conductors  embedded  in  a  multilayered  medium  above  a  ground  plane. 
Let  Qt  denote  the  region  exterior  to  the  conductors.  For  the  finite  mathematics  techniques 
[17]  and  [29],  the  unbounded  outer  region  must  be  truncated  and  enclosed  with  an 
outer  boundary  Tz  In  Chapter  2,  we  rederived  the  m1*1  order  BGT  operator  for  Laplace's 
equation.  Due  to  the  existence  of  a  ground  plane,  the  general  solution  to  Laplace's  equation 
is  slightly  different  from  (2.39).  It  reads  as 

w 

Ean 

—  cos  n<J) 

n=l  P" 


(4.1) 
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Figure  4.1.  Multiconductor  transmission  line  in  a  multilayered  dielectric 
region  above  a  ground  plane. 
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where  u  is  the  potential  function.  Hence,  B2,  the  second-order  BGT  operator,  is  also 
slightly  different  from  the  one  derived  in  Chapter  2.  For  the  asymptotic  form  of  the 
solution  to  Laplace's  Equation  (4.1),  the  B2  operator  takes  the  form 


As  will  be  seen  later,  the  finite  element  formulation  of  the  problem  involves  an  integral  over 
the  outer  boundary,  T2,  and  the  normal  derivative  of  u  appears  in  this  integrand.  Hence, 

for  our  purposes,  it  is  more  desirable  to  find  an  asymptotic  representation  for  the  normal 
derivative  of  u  rather  than  make  direct  use  of  the  operator  B2.  For  a  circular  outer 
boundary,  the  normal  derivative  will  simply  be  the  radial  one.  Using  (4.2)  along  with 
Laplace's  equation  to  exchange  the  second-order  derivative  in  p,  upp,  with  the  second- 
order  angular  derivative,  u<j)<j),  we  obtain  the  desired  asymptotic  boundary  condition: 

up  =  a(p)  u  +  P(p)  u^  (4.3) 


where  a(p)  and  p(p)  are  given  by 

«(P>  =  M.4) 

PCp)  -  ^  (4.5) 


For  a  general,  noncircular  outer  boundary  it  is  necessary  to  generalize  (4.3)  and  derive  an 
expression  for  the  normal  derivative  operator  un  in  the  local  coordinate  system  (t,n)  where  t 
and  n  are  tangent  and  normal  to  the  boundary,  respectively.  An  approximate  expression  for 
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this  derivative  can  be  obtained  by  following  the  procedure  described  in  Chapter  3  and  is 
given  by 

|“  =  cx  u+y  ut+p  utt  (4.6) 


where 


a  =  (x0  sin0o  -  y0  cos0o) - - 

3p2 


(4.7) 


Y  =' 


1  2  2 
-t  (x0  sin0o  -y0  cos0o)  +-sin20o  (yQ  -x0)  +x0y0  cos20o 


(4.8) 


-Y-1 


P  =  sin0o  -y0  cos0o)  ^  2 


(4.9) 


and  0o,  xo,  yo.  and  t  are  as  shown  in  Figure  3.2. 

4.3  Finite  Element  Implementation  of  the  Asymptotic  Boundary  Condition 

The  problem  at  hand  is  to  solve  for  the  potential  u  satisfying  the  Laplace  equation: 

V(eVu)  =  0  (4.10) 

Multiplying  (4.10)  by  a  testing  function  v  and  integrating  over  the  domain  of  the  problem 
Qj,  we  obtain 
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j  v  V-(eVu)  ds  =  0 
Jnr 


(4.11) 


From  Green's  identity  we  have 


J  v  V-(eVu)  ds  =  -J  eVu-Vv  ds  +  J  v  e^ 
^r  r2 


dl 


(4.12) 


Inserting  the  above  in  (4.1 1),  we  obtain 


3u 


J*  eVu-Vv  ds  =  J  v  e^j-  dl 

t*. 


(4.13) 


In  the  finite  element  formulation,  one  sets  up  a  mesh  in  the  region  typically  using 
triangular  elements.  The  edges  of  the  outermost  elements  prescribe  T2-  Hence, 
considering  one  element  at  a  time,  the  asymptotic  boundary  condition  given  in  (4.6)  may  be 
incorporated  into  (4.13)  to  yield 


eVu-Vv  ds  = 


l 


v  e  ( a  u  +  y  ut  +  {3  utt)  dl 


(4.14) 


Since  e  is  constant  over  each  element,  we  can  integrate  (4.14)  by  parts  to  obtain 


eVu-Vv  ds  =  e  (a  vu  +  y  vu  -  B  vtut  -  £  vut)  dl 
Jar  Jr2 


(4.15) 


where 
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—  3  (  2  \ 

C  =  (xosin0o-yocos9o)  (t  +  xocos0o+yosin0o) - - 

V  3  p  ) 


(4.16) 


The  form  given  in  (4.15)  is  well-suited  for  numerical  implementation  for  any  arbitrary 
conductor  configuration  enclosed  by  any  arbitrarily- shaped  outer  boundary. 


4.4  Numerical  Results 
4.4.1  One  conductor 


The  microstrip  line,  shown  in  Figure  4.2,  is  enclosed  by  a  rectangular  outer 
boundary.  We  first  solved  the  potential  problem  by  applying  the  asymptotic  boundary 
condition,  given  in  Equation  (4.6),  on  a  rectangular  outer  boundary.  Next,  we  introduced 
a  p.e.c.  shield  at  the  outer  boundary  and  solved  the  problem  once  again  using  the  same 
mesh.  As  Table  4.1  indicates  the  relative  error  between  the  asymptotic  boundary  condition 
and  the  published  results  [45,  52,  53]  is  between  0.053  and  2.56  percent,  whereas  the 
error  between  the  shield  and  the  published  results  is  between  14.55  and  30.04  percent.  It 
is  worth  mentioning  that  the  distance  d  from  the  microstrip  line  to  the  outer  boundary  in  the 
x-  and  y-directions  was  chosen  to  be  the  same  just  for  convenience,  but  in  principle  it  needs 
not  to  be  the  same. 


4.4.2  Two  conductors 


Two  coupled  microstrips,  shown  in  Figure  4.3,  are  enclosed  by  a  rectangular  outer 
boundary.  Table  4.2  presents  some  results  for  the  same  problem  that  have  been  published 
elsewhere  [55],  together  with  those  obtained  by  using  a  p.e.c.  shield  and  the  asymptotic 
boundary  condition  in  (4.6).  The  relative  error  in  the  self-terms  is  0.27  percent  for  the 
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Table  4.1a.  Characteristic  impedance  in  Ohms  for  the  microstrip  line  of  Figure  4.2 
(er  =  6.0). 


w/h 

ABC 

Shield 

Ref.  [45] 

Ref.  [52] 

Ref.  [53] 

Bound. 

location 

Error 

ABC- 

[53] 

Error 

Shield- 

[53] 

0.4 

90.797 

62.892 

92.278 

91.172 

89.909 

0.65 

0.987% 

30.049% 

0.7 

73.000 

53.543 

73.962 

73.613 

71.995 

0.75 

1.396% 

25.630% 

1.0 

62.531 

45.650 

62.811 

62.713 

60.970 

0.75 

2.560% 

25.127% 

2.0 

41.922 

34.778 

42.998 

43.149 

41.510 

1.20 

0.992% 

16.218% 

4.0 

26.047 

22.018 

26.971 

27.301 

26.027 

1.50 

0.077% 

15.403% 

Table  4.1b.  Characteristic  impedance  in  Ohms  for  the  microstrip  line  of  Figure  4.2 

(er  =  9.5). 


w/h 

ABC 

Shield 

Ref.[45] 

Ref.  [5  2] 

Ref.  [5  3] 

Bound. 

location 

Error 

ABC- 

[53] 

Error 

Shield- 

[53] 

0.4 

73.380 

51.513 

74.897 

73.702 

73.290 

0.65 

0.123% 

29.713% 

0.7 

58.955 

43.841 

59.910 

59.379 

58.502 

0.75 

0.774% 

25.061% 

1.0 

50.453 

37.395 

50.810 

50.501 

49.431 

0.75 

2.067% 

24.350% 

2.0 

33.766 

28.322 

34.674 

34.592 

33.493 

1.20 

0.815% 

15.439% 

4.0 

20.917 

17.863 

21.668 

21.763 

20.906 

1.50 

0.053% 

14.556% 

I 

I 

I 

I 


Outer  Boundary 

(PEC  or  ABC) 


Figure  4.2.  Microstrip  line- 
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asymptotic  boundary  condition  and  18.31  percent  for  the  shield,  and  the  error  in  the  mutual 
terms  is  5.2  percent  for  the  asymptotic  boundary  condition  and  44.59  percent  for  the 
shield. 


Table  4.2.  Capacitance  matrix  for  the  coupled  microstrips  of  Figure  4.3. 


C(i,j) 

ABC 

Shield 

Reference  [55] 

Error 

ABC-[55] 

Error 

Shield- 

155] 

C(l,l) 

0.9249  x  10-10 

0.1091  x  10-0 

0.9224  x  10-1° 

0.271% 

18.310% 

C(l,2) 

-0.8061  x  10-1 1 

-0.4712  x  10-H 

-0.8504  x  10-H 

5.203% 

44.591% 

C(2,l) 

-0.8061  x  10-1 1 

-0.4712  x  10-H 

-0.8504  x  10-H 

5.203% 

44.591% 

C(2,2) 

0.9249  x  10-10 

0.1091  x  10*9 

0.9224  x  10-10 

0.271% 

18.310% 

4.4.3  Six  conductors 

The  six  conductor  system,  shown  in  Figure  4.4,  is  enclosed,  once  again,  by  a 
rectangular  outer  boundary.  To  the  best  of  our  knowledge,  there  are  no  published  results 
for  this  configuration.  However,  we  have  compared  our  results  with  those  derived  by 
using  the  computer  program  developed  by  Harms  et  al.  [56],  which  uses  an  integral 
equation  formulation  and  an  iterative  method  of  solution.  As  Table  4.3  indicates,  the 
relative  error  for  the  capacitance  matrix  is  between  0.84  and  14.52  percent  for  the 
asymptotic  boundary  condition  and  between  3.96  and  74.10  percent  for  the  p.e.c  shield. 
That  the  distance  D  (Figure  4.4)  from  the  microstrip  line  to  the  outer  boundary  in  the  x-  and 
y-directions  was  chosen  to  be  the  same  just  for  convenience,  but  in  principle  it  needs  not  to 
be  the  same. 
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Table  4.3.  Capacitance  matrix  for  the  six-condutor  structure  of  Figure  4.4. 


C(i,j) 

Iterative  [56] 

ABC 

Shield 

Error 

ABC- 

[56] 

Error 

Shield- [5  6] 

C(l,l) 

0.668  x  10- ^ 

0.686  x  10- 

0.848  x  10- ■ 10 

2.620% 

26.83% 

C(l,2) 

-0.279  x  10- ■ lu 

-0.315  x  10- 

-0.264  x  10- 1U 

13.05% 

5.340% 

C(l,3) 

-0.549x  10- 11 

-0.600  x  10- 11 

-0.371  x  10-“ 

9.240% 

32.56% 

C(l,4) 

-0.208x  10- 11 

-0.225  x  10-“ 

-0.117x10-“ 

8.320% 

43.71% 

C(l,5) 

-0.999  x  10-1^ 

-0.101  x  10“* 

-0.456  x  10- ■ “ 

0.840% 

54730% 

T(W 

-0.704  x  10'^ 

4X6027T(FT2“ 

-0.182  x  10'“ 

14.52% 

74.10% 

C(2,l) 

-0.279  x  10- iU 

-0.315  x  10*iU 

-0.264  x  10- iU 

13.05% 

5.34% 

~&20) 

0.789  x  10- 

0.848  x  10‘ltj 

0.876  x  10- 1U 

7.480% 

11.03% 

T(23T 

-0.256  x  10- 

-0.284  x  10-lti 

-0.266  x  10- iU 

11.10% 

3.960% 

C(2,4) 

-0.465  x  10- 11 

-0.487  x  10- ll 

-0.385  x  10*“ 

4.740% 

17.23%  . 

C(2,5) 

-0.173  x  10- 11 

-0.181  x  10“  f 

-0.122  x  10-“ 

4.610% 

29.03% 

W 

-0.999  x  10' ^ 

-0.101  x  io-1  1 

-0.456  x  10- 12 

0.920% 

54.30% 

"TOT 

-0.549  x  10- 11 

-0.600x  10-“ 

-0.371  x  10'“ 

9.240% 

32.56% 

TOO) 

-0.256  x  10- 10 

-0.284  x  10- W 

-0.266  x  10- 1U 

11.10% 

3.960% 

C(3,3) 

0.794  x  10- 1U 

0.855  x  10-^ 

0.874x  10-^ 

7.680% 

10.14% 

C(3,4) 

-0.254  x  10-*0 

-0.282  x  10* 10 

-0.267  x  10* 1U 

10.91% 

4.870% 

C(3,5) 

-0.465x  10- 11 

-0.487  x  10-H 

-0.385  x  10-11 

4.750% 

17.23% 

C(3,6) 

-0.208  x  10- ^ 

-0.226  x  10- 

-0.117  x  10- 

8.410% 

43.71% 

C(4,l) 

-0.208  x  10- 11 

-0.225  x  10'“ 

-0.117  x  10-1 1 

8.320% 

43.71% 

C(4,2) 

-0.465  x  10- 11 

-0.487  x  10'“ 

-0.385  x  10-“ 

4.740% 

17.23% 

C(4.3) 

-0.254  x  10- 1U 

-0.282  x  10* 

-0.267  x  10-1° 

10.91% 

4.870% 

0.794  x  10- 10 

0.855  x  10- w 

0.874  x  10-in 

7.680% 

10.14% 

C(4,5) 

-0.256  x  10- 1U 

-0.284  x  10- ^ 

-0.266  x  10-1° 

11.10% 

3.960% 

C(4,6) 

-0.549  x  10- 11 

-0.601  x  10'“ 

-0.371  x  10-“ 

9.310% 

32.56% 

0(5,1) 

-0.999  x  10- ^ 

-0.101  x  10'“ 

-0.456  x  10- “ 

0.840% 

54.30% 

C(5,2) 

-0.173  x  10'^ 

-0.181x10-“ 

-0.122  x  10-“ 

4.610% 

29.03% 

T(53y 

-0.465x  10'“ 

-0.487  x  10-H 

-0.385  x  10-“ 

4.750% 

17.23% 

C(5,4) 

-0.256  x  10"  ^ 

-0.284  x  10- iU 

-0.266  x  10-1“ 

1 1.10% 

3.960% 

C(5,5) 

0.789  x  10-i(J 

0.848  x  10“u 

0.876  x  10-1U 

7.480% 

11.03% 

0(5,6) 

-0.279  x  10*^ 

-0.316  x  10'1U 

-0.264  x  10*  10 

13.08% 

5.340% 

C(6,l) 

-0.704  x  10’“ 

-0.602  x  10-12 

-0.182  x  10-“ 

14.52% 

74.10% 

TtfaT 

-0.999  x  10- 12 

-0.101  x  10-“ 

-0.456  x  10'“ 

0.920% 

54.30% 

C(6,3) 

-0.208  x  10"1  i 

-0.226  x  10-“ 

-0.117x10-“ 

8.4 16% 

43.71% 

0(6,4) 

-0.549  x  10-^ 

-0.601  x  10- 11 

-0.371x10-“ 

9.310% 

32.56% 

C(6,5) 

-0.279  x  IFTU 

-0.316  x  10“u 

-0.264  x  10-m 

13.08% 

5.340% 

0(6,6) 

0.668  x  10* 1U 

0.685  x  10- 1U 

0.848  x  10-m 

2.550% 

26.83% 
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In  all  of  the  three  numerical  examples  considered  thus  far,  we  have  chosen  a 
rectangular  outer  boundary  because  it  facilitates  the  meshing  procedure  and  because  it  is 
conformal  to  the  structures  considered.  This  meshing  is  done  in  a  manner  such  that  none 
of  the  triangular  elements  have  more  than  one  edge  on  the  outer  boundary.  Since  the  finite 
element  scheme  considers  oi.e  element  at  a  time,  the  problem  of  the  undefined  normal  at  the 
rectangular  comers  is  circumvented  when  the  procedure  described  above  is  followed. 

In  order  to  illustrate  the  fact  that  the  formulation  described  in  this  work  is  also 
applicable  to  an  arbitrary  outer  boundary,  we  reconsider  the  two-conductor  example 
described  earlier.  Figure  4.5  illustrates  the  coupled  microstrips  problem  with  an  arbitrary 
outer  boundary.  Table  4.4  shows  that  the  results  obtained  for  this  case  are  almost  identical 
to  those  derived  with  the  rectangular  boundary. 


Table  4.4.  Capacitance  matrix  for  the  coupled  microstrips  of  Figures  4.3  and  4.5. 


ABC  with  rectangular 

ABC  with  arbitrary 

C(i,j) 

outer  boundary 

Reference  [55] 

outer  boundary 

C(l,l) 

0.9249  x  lO-™ 

0.9224  x  10- 10 

0.9284  x  10- 10 

C(l,2) 

-0.8061  x  10- 11 

-0.8504  x  10'11 

-0.8036  x  10- 11 

C(2,l) 

-0.8061  x  10- H 

-0.8504  x  10- 11 

-0.8036  x  10'11 

C(2,2) 

0.9249  x  lO-iO 

0.9224  x  lO'iO 

0.9284  x  10- 10 

In  order  to  illustrate  the  fact  that  the  asymptotic  boundary  condition  yields  a 
significant  saving  in  computer  time  and  storage  over  the  p.e.c  shield,  we  reconsider  the 
examples  of  one-  and  two-conductor  described  earlier.  Table  4.5  shows  the  characteristic 
impedance  for  the  microstrip  line  of  Figure  4.2  using  a  p.e.c.  shield  as  a  boundary 
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condition  for  different  values  of  the  boundary  location  (d).  As  Table  4.5  indicates,  the 
p.e.c.  shield  has  to  be  placed  at  a  distance  d=4,  instead  of  d=l  as  for  the  asymptotic 
boundary  condition,  in  order  to  achieve  an  acceptable  accuracy.  Using  a  density  of  64 
nodes/unit  square,  the  configurations  of  Figure  4.2  required  a  mesh  size  of  384  and  2880 
nodes  if  an  asymptotic  boundary  condition  and  a  p.e.c.  shield  were,  respectively,  used. 
For  a  sparse  matrix  solver  the  factorization  time  is  proportional  to  N2  while  the  matrix  fill 
time  is  proportional  to  N,  where  N  is  the  total  number  of  nodes.  Hence,  the  use  of  the 
asymptotic  boundary  condition  instead  of  a  p.e.c.  shield  for  the  one-conductor 
configuration  of  Figure  4.2  resulted  in  a  reduction  factor  of  56.25  in  the  factorization  time 
and  7.5  in  the  matrix  fill  time.  In  addition,  the  storage  requirement  is  also  reduced  by  at 
least  a  factor  of  7.5. 

Table  4.5.  Characteristic  impedance  in  Ohms  for  the  microstrip  line  of  Figure  4.2  using  a 
p.e.c.  shield  as  an  outer  boundary  (er  =  6.0  and  w/h  =  1.0). 


boundary 
location  (d) 

Shield 

Reference  [53] 

Error 

Shield- [53] 

0.75 

45.650 

60.970 

25.12% 

3.0 

57.33 

5.96% 

4.0 

58.869 

3.44% 

Figure  4.6  shows  the  coupled  microstrips  with  the  outer  boundary  placed  at  a 
distance  equal  to  7.5,  instead  of  1.5  as  in  Figure  4.3,  in  the  x-d’rection  and  8.0,  instead  of 
1.7  as  in  Figure  4.3,  in  the  y-direction.  As  Table  4.6  and  Figure  4.6  indicate,  the  coupled 
microstrips  configuration  required,  using  a  density  of  64  nodes/unit  square,  a  mesh  size  of 
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2120  and  12305  nodes  if  an  asymptotic  boundary  condition  and  a  p.e.c.  shield  were, 
respectively,  used.  For  a  sparse  matrix  solver  the  factorization  time  is  proportional  to  N2 
while  the  matrix  fill  time  is  proportional  to  N,  where  N  is  the  total  number  of  nodes. 
Hence,  the  use  of  the  asymptotic  boundary  condition  instead  of  a  p.e.c.  shield  for  the 
coupled  microstrips  configuration  of  Figure  4.3  and  4.6  resulted  in  a  reduction  factor  of 
33.68  in  the  factorization  time  and  5.8  in  the  matrix  fill  time.  In  addition,  the  storage 
requirement  is  also  reduced  by  at  least  a  factor  of  5.8.  Obviously,  there  is  great  advantage 
in  using  an  asymptotic  boundary  condition. 

Table  4.6.  Capacitance  matrix  for  the  coupled  microstrips  of  Figure  4.6. 


C(i,j) 

Shield 

Reference  [55] 

Error 

Shield-[55] 

C(l,l) 

o 

• 

o 

rH 

X 

m 

00 

00 

o 

0.9224  x  10-10 

4.1% 

C(l,2) 

-0.7982  x  10- 11 

-0.8504  x  10- 11 

6.14% 

C(2,l) 

-0.7982  x  10*11 

-0.8504  x  10-H 

6.14% 

C(2,2) 

0.8845  x  10- 10 

0.9224  x  10-10 

4.1% 

Finally,  we  point  out  that  no  special  treatment  is  needed  at  the  dielectric  interfaces 
because  in  the  finite  element  formulation  the  medium  is  modeled  as  being  homogeneous 
within  each  element  and,  consequently,  the  line  integral  in  Equation  (4.15)  is  always 
confined  to  within  a  homogeneous  region  inside  the  element. 

The  term  "Error"  that  appeared  several  times  in  the  Tables  should  not  be  interpreted 
as  an  absolute  error  but  rather  as  a  difference  between  our  results  and  the  published  ones. 
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Figure  4.5.  Coupled  microstrips  with  arbitrary  outer  boundary. 


Figure  4.6.  Coupled  microstrips  with  rectangular  outer  boundary. 
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4.5  Conclusions 

In  this  chapter,  we  have  shown  how  an  asymptotic  boundary  condition  for  quasi¬ 
static  fields  can  be  applied  to  a  potential  field  at  distances  quite  close  to  a  transmission  line 
configuration  to  derive  an  FEM-based  solution  in  a  numerically  efficient  manner.  It  is 
evident  from  the  numerical  results  that  the  asymptotic  boundary  condition  consistently 
yields  more  accurate  results  than  those  obtainable  with  a  perfectly  conducting  shield  placed 
at  the  same  location.  However,  in  some  situations,  the  accuracy  obtained  with  the 
approximate  ABC  presented  in  this  chapter  may  still  not  be  adequate,  as  for  instance  in  the 
off-diagonal  terms  of  the  capacitance  matrix  in  the  six-conductor  configuration.  This  aspect 
of  the  problem  will  be  investigated  further  in  Chapter  6  where  an  improved  version  of  the 
asymptotic  boundary  condition  will  be  derived.  In  the  next  chapter,  the  three-dimensional 
version  of  the  asymptotic  boundary  condition  will  be  addressed. 
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CHAPTER  S 

ASYMPTOTIC  BOUNDARY  CONDITIONS  FOR  QUASI-TEM 
ANALYSIS  OF  THREE-DIMENSIONAL  TRANSMISSION  LINE 

DISCONTINUITIES 


5.1  Introduction 

In  Chapter  4,  an  asymptotic  boundary  condition  was  introduced  and  applied  to  two- 
dimensional  transmission  line  structures.  Typically,  a  printed  circuit  board  contains  not 
only  uniform  transmission  line  etches  that  are  essentially  invariant  in  the  longitudinal 
direction,  but  also  chip  sockets  and  connectors  for  interboard  communication  that  can  not 
be  modeled  as  uniform  lines.  Furthermore,  the  transmission  lines  themselves  may  have 
various  discontinuities  such  as  bends,  changes  in  width,  open  circuits,  gaps  and  steps.  In 
recent  years,  there  has  been  an  increasing  interest  in  modeling  such  discontinuities,  and  a 
number  of  papers  [4]  and  [57]-[72]  have  been  written  on  this  subject.  In  most  of  these 
papers,  the  integral  equation  technique  has  been  used  to  study  planar  conductors  and 
structures  containing  a  homogeneous  dielectric  with  planar  interfaces.  Castillo  [4]  has  used 
the  Finite  element  method  (FEM),  which  can  handle  any  arbitrary  configuration  of 
conductors  and  dielectrics.  When  using  the  FEM,  one  needs  to  deal  with  the  practical 
problem  of  mesh  truncation  and  the  large  number  of  mesh  nodes.  Similar  to  the  two- 
dimensional  problems,  the  most  widely-used  approach  for  dealing  with  the  mesh  truncation 
problem  for  the  three-dimensional  geometry  is  to  place  a  fictitious,  box-type  conducting 
enclosure  sufficiently  far  from  the  structure  [4],  This  approach,  which  assumes  that  the 
field  decays  significantly  before  reaching  the  outer  boundary,  typically  results  in  an 
undesirably  large  mesh,  especially  for  three-dimensional  geometries.  In  Chapter  4,  we 


74 


have  introduced  an  asymptotic  boundary  condition  (ABC),  which  provided  us  with  an 
efficient  means  for  dealing  with  open  region  two-dimensional  microwave  transmission  line 
problems  in  the  quasi-static  regime.  The  usefulness  of  the  ABC  for  obtaining  an  accurate 
solution  to  a  problem  with  a  reasonable  number  of  node  points  was  demonstrated  in  that 
chapter.  For  three-dimensional  problems,  where  the  total  number  of  mesh  points  is  usually 
large,  it  is  expected  that  the  availability  of  an  accurate  ABC  will  play  an  even  more  crucial 
role  in  the  solution  of  practical  problems. 

In  this  chapter,  we  use  a  similar,  but  more  elaborate  approach  than  in  Chapter  4  to 
derive  an  asymptotic  boundary  condition  for  three-dimensional  open  region  problems  in  the 
quasi-static  regime.  Once  again,  this  asymptotic  boundary  condition  enables  us  to  bring  the 
outer  boundary  much  closer  to  the  structure  than  would  be  possible  with  the  p.e.c.  artificial 
boundary.  In  order  to  reduce  the  number  of  unknowns  as  much  as  possible,  we  have 
chosen  an  outer  boundary  in  the  shape  of  a  parallelepiped  (we  will  refer  to  it  in  this  chapter 
as  a  box  for  the  sake  of  brevity)  because  it  is  the  most  conformable  to  the  structures 
considered. 

5.2  Derivation  of  the  Three-Dimensional  Asymptotic  Boundary  Conditions 

The  asymptotic  or  the  absorbing  boundary  condition  has  seen  an  increasing  use  in 
connection  with  the  partial  differential  equation  (PDE)  techniques  for  solving  open  region 
electromagnetic  problems  because  it  preserves  the  sparsity  of  the  discretized  PDE  matrix 
[16]-[29]  and  [32].  In  this  section,  we  use  a  similar  approach  to  the  one  in  Chapter  4  to 
derive  an  asymptotic  boundary  condition  for  three-dimensional  quasi-static  problems. 

Consider  the  three-dimensional  open  region  problem  consisting  of  an  arbitrarily- 
shaped  discontinuity  embedded  in  a  multilayered  medium  above  a  ground  plane  shown  in 
Figure  5.1.  Let  Qq-  be  the  region  exterior  to  the  conductors  and  I~2  be  the  outer  boundary. 


the  pr< 
.  a  mi 
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Our  objective  is  to  derive  an  operator  which,  when  applied  on  the  outer  boundary,  makes 
the  field  emulate  the  asymptotic  behavior  at  infinity,  and  thus  yields  an  accurate  result  for 
the  interior  region  with  only  a  moderate  number  of  nodes.  Equivalendy,  we  accomplish  this 
task  by  imposing  an  asymptotic  boundary  condition  (ABC)  on  the  field  on  the  outer 


boundary. 


The  boundary  value  problem  to  be  solved  can  be  expressed  by  the  set  of  equations: 


V-(eVu)  =0  in 


u  =  gi  on  the  im  conductor 


Bmu  =  0  on  r2 


where  u  is  the  electrostatic  potential,  gj  is  the  potential  on  the  conductors,  and  Bm  is  the 
mth  order  asymptotic  boundary  operator. 

For  large  r,  the  general  solution  of  the  Laplace  equation  in  spherical  coordinates  can 
be  written  in  terms  of  spherical  harmonics  and  powers  of  r  as  [62] 


u(r.0.<t>)  =  X  X  % 

1=0  m=-l  r 


where  Yim(0,<5))  are  spherical  harmonics.  Equation  (5.4)  can  be  rewritten  in  the  following 


u(r,0,<j>)  =  X  TT  Fin»(0»*) 


(5.5) 
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where 

l 

Flm(e.4>)=SBlmYto(6’‘t,)  <5'6) 

m=-l 


or  more  explicitly  as 


u(r,e,4>)  =  J  Flm(e,<D)  +  4  F2m(e,4>)  +  -J  F3m(e,4>)  + 


(5.7) 


We  will  now  derive  a  set  of  boundary  condition  operators  that  can  be  applied  on  the 
artificial  boundary.  From  (5.7),  we  note 


(5.8) 


We  then  define  the  first-order  operator  B  i  as 


Bju 


9u 

dr 


u 

r 


(5.9) 


The  second-order  operator  can  be  obtained  by  letting  v  =  B  i  u  and  by  observing  that 


S+t-# 


(5.10) 


The  second-order  operator  B2  can  thus  be  defined  as 
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3  Y9u 
r  A  9r  + 


(5.11) 


In  general,  it  can  be  shown  that  the  m1*1  order  asymptotic  boundary  condition  operator  can 
be  written  as 


(5.12) 


The  boundary  contribution  in  the  finite  element  formulation  enters  into  a  surface  integral 
representation  over  the  outer  boundary,  T2,  where  the  integrand  is  the  product  of  a  testing 
function  and  the  normal  derivative  of  u.  As  a  consequence,  the  asymptotic  boundary 
condition  needs  to  be  imposed  on  the  normal  derivative  of  u.  For  a  spherical  outer 
boundary,  the  normal  derivative  is  simply  the  radial  one.  Using  (5.1 1)  in  conjunction  with 
Laplace's  equation  in  the  spherical  coordinates,  we  obtain  the  following  asymptotic 
boundary  condition  operator 

ur  =  a(r)u  +  p(r)u0  +  y(r)uee  +  q(r)u0<5  (5.13) 


where 


P(r)  = 


cot© 

ir 


(5.14) 

(5.15) 


(5.16) 
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$(r)  = 


1 

2r  sin20 


(5.17) 


As  indicated  earlier,  we  need  to  use  a  conformable  outer  boundary  in  order  to  minimize  the 
number  of  node  points  as  much  as  possible.  For  three-dimensional  transmission  line 
structures,  this  conformable  outer  boundary  should  be  a  box.  It  is,  therefore,  necessary  to 
derive  the  appropriate  normal  derivative  expressions  for  the  different  faces  of  the  box 
representing  the  outer  boundary. 

For  the  x=constant  face  of  the  box,  the  normal  derivative  is  plus  or  minus  ux. 
Using  the  Chain  rule,  we  can  write  ux  as 

_  3u  3r  3u  30  3u3$  , 

Ux  3r  3x  +  30  3x  +  3$  3x 


Using  the  relations  between  the  Cartesian  and  spherical  coordinates,  Equation  (5.18)  can  be 
rewritten  as 


3u 


f  zx  > 

i  f  y  \ 

l r  J  30  ( 

v  r2p  j 

I 

(5.19) 


where  p  =  (x^  +  y2)l/2  ancj  r  =  (x2  +  y2  +  z2)l/2 

Using  (5.12)  -  (5.19),  the  Chain  rule,  and  the  relations  between  the  angular  and  the 
tangential  derivatives,  we  obtain  a  final  expression  for  the  normal  derivative  on  the  face 
where  x=constant  and  ux  is 

ux  =  aj(x,y,z)u  +  Pj(x,y,z)uz  +  ^(x.y.zlu^ 


+  £j(x,y,z)uy  +  ri^x^.zluyy  +  ^(x,y,z)uyZ 


(5.20) 
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where 


a^x.y.z)  =  - 


(5.21) 


Pi(x,y,z)  =  -z 


3x3+4xy2 


(5.22) 


Yi(x,y,z) 


X3  +  XV2 


2r 


(5.23) 


^(x,y,z) 


4x3yz2+  3xy3z2-xy(p4  +  2rp2) 
2r2p4 


(5.24) 


rii(x,y,z)  = 


z2y2x  +  x3r2 

2r2p2 


(5.25) 


^i(x,y,z) 


zxy 


(5.26) 


On  the  y=constant  face,  the  normal  derivative  is  plus  or  minus  uy.  Invoking  the  Chain 
rule,  once  again,  we  can  express  uy  as 


_  9u  9r  9u  99  9u  9p 
Uy  9r  9y  +  90  9y  +  9p  9y 


(5.27) 


Following  the  same  procedure  as  before,  we  obtain  the  following  expression  for  uv 
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Uy  =  a2(x,y,z-)u  +  p2(x,y,z)uz  +  Y2^x»y.z)uz2 

+  £2(x,y,z)ux  +  Tj2(x,y,z)uxx  +  ^2(x,y,z)uxz 


(5.28) 


where 


a2(x,y,z)  =  - 


y 


(5.29) 


p2(x,y,z) 


3y3+  4yx2 
2r2p2 


(5.30) 


y2(x,y,z) 


y3+  yx 

2r2 


2 


(5.31) 


k&.y.z) 


4y3xz2+  3yx:V-xy(p4  +  2r2pz) 


3_2 


2.2\ 


2r2p4 


(5.32) 


Tl2(x,y,z)  = 


zVy.yV 

2r2p2 


(5.33) 


^2(x«y.z) 


zxy 


(5.34) 


Similarly,  for  the  z=constant  face,  we  can  obtain  an  expression  for  uz,  which  reads: 


uz  =  a3(x,y,z)u  +  (33(x,y,z)ux  +  y3(x,y,z)uxx 

+  £3(x,y,z)uy  +  Ti3(x,y,z)uyy  +  £3(x,y,z)uxy 


(5.35) 


where 
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\  z 

a3u,y,z,  =  — - 

r 


P3(x,y,z)  = 

2r 


2  3  2  2 

/  \  x  z  +  y  zr 

y3(x,y,z)  = - — - 

2rp2 


C3(x,y,z)  = 

2? 


,  N  y2z3  +  x2zr2 

Ti3(x,y,z)  = - r— - 

2r  P 


43(x,y,z)  =  — 


(5.36) 


(5.37) 


(5.38) 


(5.39) 


(5.40) 


(5.41) 


It  is  evident  that  one  needs  to  use  the  appropriate  normal  derivative  expression  on  different 
faces  of  the  box-shaped  outer  boundary.  As  can  be  seen  from  Equations  (5.12)-(5.41),  it 
is  much  easier  to  choose  a  spherical  outer  boundary  where  the  normal  derivative  is  simply 
ur.  However,  for  the  purpose  of  truncating  the  unbounded  region  surrounding  the 
transmission  lines  in  an  efficient  manner,  one  needs  to  use  a  conformable  outer  boundary 
which,  as  mentioned  earlier,  is  typically  a  box- shaped  surface  for  three-dimensional 
transmission  line  structures.  The  asymptotic  boundary  condition  expressions  that  we  have 
just  derived  will  be  implemented  in  the  finite  element  scheme  in  the  next  section. 
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5.3  Finite  Element  Implementation  of  the  Asymptotic  Boundary  Condition 


The  finite  element  method  and  its  implementation  are  well-documented  in  the 
literature  [1],  [3]-[6],  and  [30]-[31].  In  this  section,  we  will  be  concerned  with  the 
implementation  of  the  asymptotic  boundary  condition  in  the  finite  element  method.  As 
indicated  earlier,  the  region  of  interest,  Qj,  is  bounded  by  an  artificial  boundary,  f2  ,  to 
limit  the  number  of  unknowns.  Over  the  bounded  region,  the  Laplace  equation  is  solved  at 
a  finite  number  of  grid  points.  This  equation  is  discretized  through  the  use  of  a  weak  form 
of  variational  representation. 

Multiplying  the  Laplace  Equation  (5.1)  by  a  testing  function  f  and  integrating  over 
the  volume  O7,  result  in 


f  fV-(eVu)  dv  =  0  (5.42) 

Jar 


Using  the  Green's  second  identity,  we  obtain 


f  fV-(eVu)  dv  =  -  j"  eVu-Vf  dv  +  f  fe  3^-  ds  (5.43) 

Jop  JOj  J  r2  °n 


Clearly,  the  second  term  of  the  right-hand  side  of  (5.43)  is  the  boundary  integral 
contribution  which  is  usually  neglected  when  a  p.e.c.  outer  boundary  is  used  [4], 

Inserting  (5.43)  into  (5.42),  we  get 

f  eVu-Vf  dv  =  f  fe  3^  ds  (5.44) 

J  r2  ™ 


The  region  is  discretized  into  tetrahedral  elements.  The  triangular  faces  of  the 
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outermost  elements  make  up  the  outer  boundary  r2.  In  finite  element  formulation,  the 
implementation  of  (5.44)  is  carried  out  on  an  element-by-element  basis.  For  all  but  the 
outermost  tetrahedral  elements,  the  right-hand  side  of  Equation  (5.44)  is  zero.  The 
asymptotic  boundary  condition  is  needed  to  treat  those  outermost  elements. 

For  those  elements  having  a  face  on  the  surface  prescribed  by  x=constant,  where 
the  outward  normal  is  in  the  plus  or  minus  x-direction,  the  asymptotic  boundary  condition 
given  in  (5.20)  may  be  incorporated  into  (5.44)  to  yield 

I  eVu-Vfdv  =±|  fe(  a1(x,y,z)u  +  P^x.y.zluj  +  y^x.y.zlu^ 

Jar  Jr2 

+  Ci(x,y,z)uy  +  Tjjtx.y.zlUyy  +  ^1(x,y,z)uyz)  dydz  (5.45) 


Similarly,  for  the  elements  having  a  face  on  the  surface  prescribed  by  y=constant,  the 
asymptotic  boundary  condition  given  in  (5.28)  is  used  in  (5.44)  to  give 

j  eVu-Vf  dv  =  ±  I  fe(  a2(x,y,z)u  +  P2(x,y,z)uz  +  y2(x,y,z)uzz 
Jfir  Jr2 

+  ^2(x,y,z)ux  +  rt2(x,y,z)uxx  +  £2(x,y,z)uxz)  dxdz  (5.46) 


For  the  elements  having  a  face  on  the  surface  prescribed  by  z=constant,  where  the  outward 
normal  is  in  the  plus  z-direction,  (5.35)  is  incorporated  into  (5.44)  to  yield 

eVu-Vf  dv  =  fe  (a3(x,y,z)u  +  p3(x,y,z)ux  +  y3(x,y,z)uxx 

Jar  Jr2 

+  ^3(x,y,z)uy  +  ri3(x,y,z)uyy  +  ^3(x,y,z)uxy)  dxdy  (5.47) 


Therefore,  for  those  elements  which  do  not  share  a  face  with  T 2,  the  boundary  contribution 
is  zero.  However,  for  the  outermost  elements,  depending  on  the  outward  normal,  the 
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appropriate  asymptotic  boundary  condition  has  to  be  used  (i.e.,  either  the  right-hand  side  of 
(5.45),  (5.46),  or  (5.47)  has  to  be  added  to  the  corresponding  element  matrix). 

5.4  Numerical  Results 

A  rectangular  section  of  a  microstrip  transmission  line  of  length  L,  width  W,  and 
height  H  above  the  ground  plane  is  shown  in  Figure  5.2.  The  outer  boundary  F2  was 
chosen  to  have  the  shape  of  a  box.  Using  the  same  mesh,  we  have  solved  the  potential 
problem  twice,  first  by  applying  the  asymptotic  boundary  condition  on  the  outer  boundary, 
and  second  by  placing  a  perfect  electric  conducting  shield  at  the  same  location.  After 
solving  for  the  electrostatic  potential,  we  computed  the  normalized  capacitance  CH/e(area) 
for  both  cases.  Tables  5.1  and  5.2  show  the  results  of  computation  for  the  normalized 
capacitance  for  different  values  of  L/W  and  for  three  dielectric  constants  (Er  =  1 .0,  6.0, 
9.6).  As  Tables  5.1  and  5.2  indicate,  the  asymptotic  boundary  condition  yields  more 
accurate  results  than  those  obtainable  with  a  perfectly  conducting  shield  [60].  Clearly,  for 
this  problem  there  is  a  distinct  advantage  in  using  an  asymptotic  boundary  condition  in 
place  of  a  p.e.c.  shield. 

(  CH  'n  l  pj 

Table  5.1.  Normalized  capacitance  (  £"^rea  j  for  =  0.2,  -^-=0.2,  Dx=Dy=Dz  =0.5. 


P.E.C. 

ABC 

er 

Shield 

(Present  Method) 

Reference  [60] 

1.0 

1.34 

3.73 

3.5 

6.0 

1.04 

2.25 

2.2 

9.6 

1.02 

2.12 

2.1 
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r  CH  \ 


H 


Table  5.2.  Normalized  capacitance  ^  £  J  for  =  1.0,  -^=1.0,  Dx=D>=r»z  =2.0 


P.E.C. 

ABC 

Er 

Shield 

(Present  Method) 

Reference  [60] 

1.0 

2.19 

4.90 

5.0 

6.0 

1.36 

3.11 

3.4 

9.6 

1.31 

2.84 

2.9 

Finally,  we  would  like  to  mention  that  the  boundary  surface  integral  is  always 
confined  within  a  homogeneous  region  inside  the  element  because  the  finite  element  method 
models  the  medium  to  be  homogeneous  within  each  element;  consequently,  no  special 
treatment  is  needed  at  the  dielectric  interfaces. 

5.5  Conclusions 

Starting  from  the  general  solution  of  the  Laplace  equation  in  spherical  coordinates, 
we  derived  a  set  of  asymptotic  boundary  conditions  for  three-dimensional  quasi-static 
problems  for  a  spherical  outer  boundary.  The  second-order  boundary  condition  was  then 
generalized  to  a  box-shaped  outer  boundary  and  implemented  in  the  finite  element  method 
to  solve  the  potential  problem  of  a  rectangular  microstrip  patch.  The  numerical  results 
show  that  the  asymptotic  boundary  conditions  yield  more  accurate  results  than  those 
obtainable  with  a  perfectly  conducting  shield  placed  at  the  same  location. 
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CHAPTER  6 

IMPROVEMENT  ON  THE  BOUNDARY  CONDITIONS: 

THE  HIGHER-ORDER  BOUNDARY  CONDITIONS 

6.1  Introduction 

The  finite  element  method  (FEM)  is  very  appealing  for  solving  open  region 
problems  due  to  its  simplicity  in  modeling  complex-shaped  structures  and  inhomogeneous 
dielectric  scatterers.  The  local,  or  the  absorbing  boundary  condition  (ABC),  makes  the 
FEM  even  more  powerful  because  it  preserves  the  sparsity  of  the  discretized  matrix. 
Furthermore,  the  absorbing  boundary  condition  operator  mimics,  to  a  certain  degree,  the 
asymptotic  behavior  of  the  wave  function  at  infinity  and  yields  reasonably  accurate  results 
in  the  interior  region  without  the  need  of  an  exorbitantly  large  number  of  mesh  points. 
Among  the  available  absorbing  boundary  conditions  [8]-[32],  the  Bayliss,  Gunzburger, 
and  Turkel  (BGT)  boundary  conditions  are  the  most  commonly  used.  In  Chapter  3,  we 
have  presented  a  generalization  of  the  original  BGT  absorbing  boundary  condition  so  as  to 
make  it  applicable  to  an  arbitrary,  rather  than  circular,  outer  boundary.  The  use  of  the 
generalized  version  of  the  BGT  enables  one  to  reduce  the  number  of  node  points 
significantly  and  to  solve  larger  sized  problems  than  had  been  possible  in  the  past.  In 
Chapter  4,  we  derived  the  static  version  of  the  BGT  boundary  condition.  As  the  numerical 
results  in  Chapter  4  show,  the  asymptotic  boundary  condition  consistently  yielded  more 
accurate  results  than  those  obtainable  with  a  perfectly  conducting  shield  placed  at  the  same 
location.  Nevertheless,  most  forms  of  the  absorbing  boundary  condition  operators  have 
been  based  on  the  use  of  only  the  first  few  terms  of  the  asymptotic  representation  of  the 
solution  to  the  differential  equation.  Typically,  the  first  two  terms  of  the  series  are  used  to 
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obtain  the  desired  absorbing  boundary  condition  operator.  In  [24],  [27],  and  [28],  it  was 
demonstrated  that  while  the  absorbing  boundary  condition,  which  is  based  on  the  first 
terms  of  the  series,  works  quite  well  for  the  lower-order  harmonics,  it  exhibits  a  significant 
error  for  the  higher-order  harmonics,  which,  in  turn,  may  cause  a  noticeable  error  in  the 
finite  element  solution. 

In  this  chapter,  we  derive  a  higher-order  absorbing  condition  which  unlike  many  of 
the  available  ABCs  that  take  into  account  only  the  lower-order  harmonics,  considers  both 
lower-  and  higher-order  harmonics.  In  common  with  the  other  available  ABCs,  the 
derivation  of  this  new  absorbing  condition  is  based  on  the  same  principle  of  the  asymptotic 
representation  of  the  solution  to  the  differential  equation.  However,  unlike  the  available 
ABCs  which  assume  that  in  the  far  region  the  solution  can  adequately  be  represented  by  the 
first  few  terms  of  the  series,  the  higher-order  absorbing  condition  requires  that  the 
asymptotic  representation  be  a  combination  of  lower-  and  higher-order  harmonics.  As  will 
be  demonstrated  later,  the  use  of  the  higher-order  ABC  results  in  a  significant  improvement 
in  the  finite  element  solution  for  a  variety  of  scatterers  and  transmission  line  structures. 

6.2  Derivation  of  the  Higher-Order  Absorbing  Boundary  Condition 

Consider  the  geometry  of  the  scattering  problem  described  in  Figure  2.3.  The 
region  £2  is  bounded  with  the  contour  T i  and  its  exterior  region  is  denoted  by  Oj.  Let  T2 
be  the  outer  contour,  which  truncates  the  open  region  Qt.  where  the  absorbing  boundary 
condition  will  be  applied.  Obviously,  one  can  now  use  a  partial  differential  equation  (PDE) 
technique  to  solve  the  approximate  problem  provided  that  the  region  bounded  by  Ti  is  not 
so  large  as  to  require  an  unmanageable  number  of  unknowns.  The  approximate  problem, 
shown  in  Figure  2.3,  is  equivalent  to 
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(6.1) 

(6.2) 

Bmu  =  0  on  r2  (6.3) 


V2u  +  k2u  =  0  in  Qt 


du 

^=«°nr, 


where  u  is  the  scattered  field,  g  is  the  contribution  from  the  incident  field,  and  Bm  is  the 
m^1  order  absorbing  boundary  condition  operator.  Unlike  the  asymptotic  representation  of 
the  scattered  field  given  in  Chapter  2,  where  only  the  first  terms  of  the  series  are  kept,  we 
suggest  the  following  asymptotic  form 


»-jkp 


Vp 


arnfo) 


V  P 


nl 


aTl2(<?>) 


n2 


.n3 


nl  <  n2  <  n3 


(6.4) 


Note  that  the  asymptotic  representation  of  the  scattered  field  given  by  (2.22)  can  be 
obtained  from  (6.4)  by  choosing  the  set  {ni,  n2,  n3)  to  be  equal  to  {0,  1,  2). 

From  (6.4),  we  can  see  that 


Taking  a  closer  look  at  the  right-hand  side  of  (6.5),  we  note  that 


du  ,u  (l  +  "' 

3F+jku  +  _ F” 


u  =  0| 


1  A 


n2+3/2 


(6.6) 


Thus  we  define  the  boundary  operator  on  the  left  side  of  (6.6)  to  be 
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(6.7) 


Note  that  if  ni  is  chosen  to  be  equal  to  1,  the  operator  Bi  is  equivalent  to  the  first-order 
BGT  operator  [29].  If  (6.7)  is  explicitly  written,  it  would  have  the  form 


.-jkp 


Bru  =  — 


.vi  -  “P  +  3T(n3  - 


nx)  +... 


(6.8) 


The  second-order  boundary  operator  can  readily  be  obtained  by  letting  v  =  B  i  u  and 
observing  that 


Note  that  if  the  set  {ni,  n2)  is  chosen  to  be  equal  to  { 1,  2),  the  operator  B2  is  equivalent 
to  the  second-order  BGT  operator  [29]. 

Similarly,  a  third-order  boundary  operator  can  be  obtained 


(6.11) 
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Note  that  if  the  set  {ni,  n2, 113}  is  chosen  to  be  equal  to  (1,  2,  3),  the  operator  B3  is 
equivalent  to  the  third-order  BGT  operator  [29]. 

As  was  indicated  in  Chapter  3,  the  boundary  integral  contribution  is  a  line  integral 
over  I"2  where  its  integrand  is  a  product  of  the  testing  function  and  the  normal  derivative  of 
the  unknown  field.  Hence,  for  our  purposes,  it  is  more  useful  to  find  an  asymptotic 
representation  for  the  normal  derivative  of  u  rather  than  make  direct  use  of  the  boundary 
condition  operators  obtained  above.  For  a  circular  outer  boundary,  the  normal  derivative 
is  simply  the  radial  one.  Using  the  B3  operator  as  given  by  (6.11)  and  the  Helmholtz 
equation  given  by  (6.1),  we  obtain  an  asymptotic  representation  for  the  radial  derivative 
that  reads 

=  a(nlf  n2,n3,p)u  +  P(n1,n2,n3,p)u0(j)  (6.12) 


where 


a(nItn2,n3,p)  =  ■“(n,n2n3+^(n1+n2+n3)+^-(n1n2+n1n3+n2n3)+-D 

+  jk^n1n2+n1n3+n2n3+2(ni+n2+n3)+-^- j-2pk2(n!+n2+n3+4)-4p2jk3)  J  / 


^nln2+nIn34-n2n3+ni+n2+n3+^  +  2pjk(n2+n2+n3+3)  -  4p2k2 


(6.13) 


and 


P(n1,n2,n3,p)  = 


— ^ni+n2+n3+ 


f)+3jk 


^n1n2+njn3+n2n3+n1+n2+n3-i~^+  2pjk(n1+n2+n3+3)  -  4p2k2 


(6.14) 
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Note  that  the  absorbing  boundary  condition  operator  shown  in  (6.12)  depends  on  the 
harmonics  (ni,  n2,  n3). 

In  deriving  (6.12),  we  have  used  the  first-order  boundary  condition  operator  Bj  to 
approximate  the  term  u<jjpp  as  follows 


The  absorbing  boundary  condition  that  we  have  derived  in  (6.12)  applies  to  a  circular  outer 
boundary  where  the  normal  derivative  of  u  is  simply  its  radial  derivative  up.  Our  next  task 

is  to  transform  the  above  absorbing  boundary  operator  into  a  form  suitable  for  arbitrary 
boundaries.  Using  the  same  transformation  and  the  same  approach  as  the  one  used  in 
Chapter  3,  we  obtain  the  following  absorbing  boundary  condition  operator  that  can  be 
applied  on  an  arbitrary  outer  boundary 


9u 

9n 


au  +  yut  +  Putt 


(6.16) 


where 
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a  =  -(xosin0o-yocos0o)£  ~^n1n2n3+^(n1+n2+n3)+-j(n1n2+n1ri3+n2n3)+-g- j 

+jk^n1n2+n1n3+n2n3+2(ni+n2+n3)+-~^-2pk2(nI+n2+n3+4)-4p2jk3)J/ 
^P^n1n2+n1n3+n2n3+n1+n2+n3+-j^  +  2p2jk(n1+n2+n3+3)  -  4p3k2j 


(6.17) 


1  0  0 

_  -t  (x0  sin0o  -yo  cos0o)  +-jsin20o  (y0  -x0)  +x0y0  cos20o 
7  - - ; - 


(6.18) 


P  =  £(xosin6o-yocos0o)(n1+n2+n3+-jj+3jkpj/[(  njn2 


J  0  0 

+n1n3+n2n3+n1+n2+n3+— )(-2jkp(n1+n2+n3+3)-4p  k 


] 


(6.19) 


and  00,  xo,  yO,  and  t  are  as  shown  in  Figure  3.2. 

The  higher-order  boundary  condition  given  by  (6.16)  is  implemented  in  the  finite 
element  scheme  in  the  same  way  as  the  modified  BGT  boundary  condition  discussed  in 
Chapter  3. 

6.3  Derivation  of  the  Higher-Order  Asymptotic  Boundary  Condition 

Consider  the  problem  of  N  arbitrarily-shaped  conductors  embedded  in  a 
multilayered  medium  above  a  ground  plane,  shown  in  Figure  4.1.  Let  Qj  denote  the 
region  exterior  to  the  conductors  and  T2  the  outer  boundary  enclosing  the  truncated  region. 
The  asymptotic  boundary  condition  should  mimic  the  asymptotic  behavior  of  the  field  at 
infinity  and  yield  reasonably  accurate  results  in  the  interior  region.  The  potential  u  must 
satisfy  Laplace's  equation  everywhere  in  £2t,  the  constant  potential  condition  on  the 
conductors,  and  the  asymptotic  boundary  condition  on  the  outer  boundary  f;. 
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Equivalently,  the  problem  can  be  described  in  terms  of  the  following  equations: 


V-(eVu)  =0  in 

(6.20) 

u  =  g;  on  the  Ith  conductor 

(6.21) 

Bmu  =  0  on  r2 

(6.22) 

where  u  is  the  electrostatic  potential  and  Bm  is  the  mfh  order  asymptotic  boundary 
condition. 

In  Chapter  4,  the  asymptotic  boundary  condition  was  based  on  the  first  two  terms 
of  the  asymptotic  representation  of  the  general  solution  to  Laplace's  equation.  In  this 
chapter,  we  suggest  the  following  asymptotic  form  for  the  potential  u 

^nl  a„2  a„3 

u  =  — ^-cosn^  +  — — cosn2<{)  +  — ^cosn3(J)  +  . . .  (6.23) 

Pn  P  P 


Note  that  the  asymptotic  representation  of  the  potential  u  given  in  Chapter  2  can  be  obtained 
from  (6.23)  by  choosing  the  set  {ni,  n2,  n3)  to  be  equal  to  { 1,  2,  3}. 

From  (6.23),  we  can  see  that 


3u  niu 


^2  £^3 

-^cosn2<J>(n1-n2)+-^j-cosn34)(n1-n3)  +  .  . . 

P  P 


(6.24) 


We,  therefore,  define  the  first-order  operator  B  l  to  be 
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5u 

®>u=^+ 


njU 

T" 


(6.25) 


Note  that  if  ni  is  chosen  to  be  equal  to  1,  the  Bi  operator  given  in  Chapter  4  can  be 
recovered  from  (6.25). 

The  second-order  boundary  operator  can  readily  be  obtained  by  letting  v  =  Bju  and 
observing  that 


dv 

dp 


+  (nJ+l£  =  of  '  ' 

P 


n3+ 2 


(6.26) 


Thus,  the  second-order  operator  is  defined  by 

b2u  ■  + + 7)“  (6-27 

Note  that  the  B2  operator  of  Chapter  4  can  be  recovered  from  (6.27)  by  choosing  the  set 
{ni,  n2)  to  be  equal  to  { 1,  2). 

The  third-order  operator  can  be  obtained  by  letting  z  =  B2v,  and  defining  B 3  to  be 
equal  to 

B3U  s  ( A  +  (n3+2)i)(^~  +  +  ni^)  (6’28 

As  explained  earlier,  we  need  to  find  an  asymptotic  representation  for  the  normal  derivative 
of  u  rather  than  make  direct  use  of  the  boundary  condition  operators  obtained  above.  For  a 
circular  outer  boundary,  the  normal  derivative  is  simply  the  radial  one.  Using  the  B3 
operator  as  given  in  (6.28)  and  making  use  of  the  Laplace  equation  to  trade  the  second- 
order  derivative  in  p,  upp,  for  the  second-order  angular  derivative,  uqp,  we  obtain  an 
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asymptotic  representation  for  the  radial  derivative  that  reads 
=  Mni.  n2,  n3,  p)u  +  £(nlf  n2,  n3,  p)u$$ 


(6.29) 


where 


ntn2n3 

A.(n^,  n2,  n3,  p)  —  / 

pui^  +  ntn3  +  n2n3) 


(6.30) 


*+■  112  4*  113 

^(n!’  °2’  °3,  ^  p(njn2  +  ntn3  +  n2n3) 


(6.31) 


In  deriving  (6.29),  we  used  the  first-order  boundary  condition  operator  Bi  to  approximate 
the  term  u^p,  that  is, 

Upoo  -  (6.32) 

The  asymptotic  boundary  condition  operator  given  by  (6.29)  is  valid  only  for  a  circular 
outer  boundary.  However,  for  the  purpose  of  reducing  the  number  of  mesh  points  as 
much  as  possible,  one  needs  to  use  a  conformable  outer  boundary.  Thus,  we  need  to 
generalize  (6.29)  and  obtain  an  operator  valid  for  an  arbitrary  outer  boundary.  Following 
the  procedure  described  in  Chapter  3,  we  obtain  an  expression  for  the  normal  derivative  un 
in  the  local  coordinate  system  (t,n)  where  t  and  n  are  tangent  and  normal  to  the  triangular 
edge  lying  on  the  outer  boundary  T2,  respectively, 

=  Wni,  n2,  n3,  p)u  +  yut  +  ^(nj,  n2,  n3,  p)u„ 


(6.33) 
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where  ut  and  utt  are  the  first-  and  second-order  tangential  derivatives,  respectively,  and 


X(nj,  n2,  n3,  p)  =  -(xosin0o-yocos9o) 


n^nj 


2 

p  (nxn2  +  n^  +  n2n3) 


(6.34) 


—  J 

5(nlt  n2,  n3,  p)  =  (xosineo-yocos0o) 


ni  -t-  n2  +  n3 


2 

p  (n^  +  n^  +  n2n3) 


(6.35) 


where  0o,  xo,  yO»  and  t  are  as  shown  in  Figure  3.2. 

6.4  Numerical  Results 

The  higher-order  boundary  condition  is  implemented  in  the  finite  element  scheme  in 
the  same  way  as  the  simple  absorbing  boundary  condition  because  they  both  yield  the  same 
form  of  absorbing  boundary  condition  operator.  To  demonstrate  the  significant 
improvement  in  the  finite  element  solution  achieved  by  using  the  higher-order  boundary 
conditions,  we  considered  both  circuits  and  scattering  problems.  For  all  the  cases 
considered,  a  conformable  outer  boundary  is  used  to  truncate  the  open  region. 

6.4.1  Digital  circuit  applications 

6.4. 1.1  Two  conductors 

Consider  the  two  coupled  microstrips  shown  in  Figure  4.3.  The  higher-order 
asymptotic  boundary  condition  operator  given  by  (6.33)  was  applied  on  a  rectangular  outer 
boundary.  Choosing  the  set  (ni,  n2,  n3 }  to  be  equal  to  {1,  2,  4),  the  finite  element 
problem  was  solved  for  the  electrostatic  potential.  Table  6.1  shows  the  capacitance  matrix 
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for  the  same  problem  that  has  been  published  elsewhere  [55],  together  with  those  obtained 
by  using  a  p.e.c.  shield,  the  asymptotic  boundary  condition  (the  method  of  Chapter  4),  and 
the  higher-order  asymptotic  boundary  condition  (present  method).  As  Table  6.1  indicates, 
while  both  the  present  method  and  the  simple  asymptotic  boundary  condition  yield  more 
accurate  results  than  those  obtainable  with  a  perfectly  conducting  shield  placed  at  the  same 
location,  the  higher-order  asymptotic  boundary  condition  results  compare  more  favorably 
with  the  published  work. 


Table  6.1.  Capacitance  matrix  for  the  coupled  microstrips  of  Figure  4.3. 


C(i,i) 

Reference  [55] 

ABC 

Higher  Order 

ABC 

Error 

ABC-[55] 

Error 

HOABC-[55] 

C  (1,1) 

0.9224  x  10- 10 

0.9249  x  10-1° 

0.9230  x  10-10 

0.271% 

0.069% 

C(l,2) 

-0.8504  x  10*11 

-0.8061  x  10*11 

-0.8377  x  10-H 

5.203% 

1.489% 

C(2,l) 

-0.8504  x  10-H 

-0.8061  x  10-H 

-0.8377  x  10-H 

5.203% 

1.489% 

C(2,2) 

0.9224  x  10-10 

0.9249  x  10-10 

0.9230  x  10-10 

0.271% 

0.069% 

6.4.1.2  Six  conductors 


Consider  the  six-conductor  system  shown  in  Figure  4.4.  As  we  indicated  in 
Chapter  4,  there  are  no  published  results  for  this  configuration.  However,  we  have 
compared  our  results  with  those  derived  by  using  the  computer  program  developed  by 
Haims  et  al.  [56],  which  uses  an  integral  equation  formulation  and  an  iterative  method  of 
solution.  For  this  configuration,  although  the  simple  ABC  of  Chapter  4  yields  more 
accurate  results  than  those  obtainable  with  a  perfectly  conducting  shield  placed  at  the  same 
location,  the  error  in  the  capacitance  matrix  is  noticeable,  especially  for  the  off-diagonal 
terms.  The  higher-order  asymptotic  boundary  condition  operator  given  by  (6.33)  was 
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applied  on  a  conformable  rectangular  outer  boundary  with  a  choice  of  the  set  {ni,  n2,  n3) 
to  be  equal  to  { 1,  5,  10}.  As  Table  6.2  indicates,  the  higher-order  asymptotic  boundary 
condition  yields  a  significant  improvement  over  the  simple  asymptotic  boundary  condition, 
especially  for  the  off-diagonal  terms  of  the  capacitance  matrix. 

Clearly,  the  improvements  brought  about  by  the  higher-order  asymptotic  boundary 
condition  are  a  direct  consequence  of  its  ability  to  incorporate  lower-  and  higher-order 
terms  through  the  choice  of  the  set  {ni,  n2,  n3).  Based  on  numerical  investigations,  it 
was  determined  that  the  optimal  choice  of  the  set  {ni,  n2,  n3)  is  { 1,  p/2,  p}  where  p  is  the 
distance  from  the  origin  to  the  middle  of  the  edge  of  the  triangular  element  residing  on  the 
outer  boundary  (see  Figure  3.2). 

6.4.2  Scattering  applications 

6.4.2.1  Perfect  electric  conductor  circular  cylinder 

Consider  a  perfect  electric  conductor  cylinder  with  a  radius  of  5\  enclosed  with  a 
circular  outer  boundary  of  radius  5.15X  shown  in  Figure  6.1.  For  a  TE  incident  wave,  the 
finite  element  method  was  used  to  obtain  an  approximate  solution  which  was  then 
compared  to  the  exact  series  solution.  Using  the  same  mesh,  the  finite  element  problem 
was  solved  twice  using  two  different  boundary  conditions:  once  by  using  the  BGT 
absorbing  boundary  condition  and  the  other  by  using  the  higher-order  absorbing  boundary 
condition.  Figure  6.2  shows  the  radar  cross-section  (RCS)  for  the  exact  series  solution, 
together  with  the  one  obtained  by  using  the  BGT  absorbing  boundary  condition  and  the 
higher-order  absorbing  boundary  condition  (present  method).  Figure  6.3  shows  the  error 
in  the  finite  element  RCS  computation.  As  Figures  6.2  and  6.3  indicate,  the  higher-order 
absorbing  boundary  condition  yields  a  significant  improvement  in  accuracy  over  the  BGT 
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Table  6.2.  Capacitance  matrix  for  the  six-conductor  structure  of  Figure  4.4. 


C(i,j) 


pgeea 


C(2.5) 


C(3,l) 


C(4,6) 


C(5,l) 


C(6,6) 


Iterative  | 

[56] 

0.668  x 

io- 1 

-0.279  x 

10-: 

-0.549  x 

10- : 

-0.208  x 

10- : 

-0.999  x 

10- 

-0.704  x 

10- 

-0.279  x 

10- 

0.789  x 

10- 1 

-0.256  x 

10- 

-0.465  x 

10- 

-0.173  x 

10- 

-0.999  x 

10- 

-0.549  x 
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Figure  6.2.  The  radar  cross-section  of  a  5  a,  cylinder  illuminated  by  a  TE 

incident  wave  and  enclosed  by  a  circular  outer  boundary  0.15  X 
away  from  the  surface  of  the  scatterer  as  shown  in  Figure  6.1. 
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Figure  6.3.  The  error  in  the  radar  cross-section  between  the  finite  element  and 
the  exact  series  solutions  for  a  5  \  cylinder  illuminated  by  a  TE 

incident  wave  and  enclosed  by  a  circular  outer  boundary  0.15  X. 
away  from  the  surface  of  the  scatterer  as  shown  in  Figure  6.1. 
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boundary  condition.  It  is  worth  mentioning  that  for  a  TM  wave  incident  on  the  same 
cylinder  the  higher-order  and  the  BGT  boundary  conditions  give  almost  the  same  result 
which  compares  very  well  with  that  for  the  exact  series  solution. 

6.4. 2.2  Perfect  electric  conductor  wedge 

Consider  a  6k  by  3 k  perfect  electric  conductor  wedge  illuminated  by  a  TM  incident 
wave  (E*  =  e*J^x  ).  In  order  to  minimize  the  number  of  mesh  points,  the  wedge  was 
enclosed  by  a  conformable  outer  boundary  1\  away  from  the  surface  of  the  scatterer 
(Figure  6.4).  The  near  field  was  calculated  on  the  outer  boundary  using  the  higher-order 
boundary  condition  (present  method),  the  BGT  boundary  condition,  and  the  method  of 
moments.  From  Figure  6.5,  we  can  deduce  that  the  use  of  the  higher-order  boundary 
condition  results  in  considerable  improvement  over  the  BGT  boundary  condition.  In  order 
to  see  more  clearly  the  difference  between  higher-order  and  the  BGT  boundary  conditions, 
we  have  plotted  the  results  obtained  by  the  two  methods  on  separate  graphs  in  Figure  6.6. 

6.4.2.3  Perfect  electric  conductor  strip 

Figure  6.7  shows  a  6k  strip  illuminated  by  a  TM  incident  wave  (E*  =  e‘j^x  ). 
Again,  for  the  purpose  of  reducing  the  number  of  mesh  points,  the  outer  boundary  was 
chosen  to  be  as  conformable  as  possible  to  the  surface  of  the  scatterer.  The  near  field  was 
calculated  on  the  outer  boundary  using  the  higher-order  boundary  condition  ,  the  BGT 
boundary  condition,  and  the  method  of  moments.  The  results  obtained  via  the  use  of  the 
higher-order  and  the  BGT  boundary  conditions  are  again  compared  to  the  method  of 
moments  results.  As  Figures  6.8  and  6.9  show,  the  higher-order  boundary  condition 
compares  more  favorably  than  the  BGT  boundary  condition  with  the  method  of  moments. 
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Figure  6.5.  Near  field  on  the  outer  boundary  for  a  6X  by  3A.  p.e.c.  wedge 

illuminated  by  a  TM  incident  wave  (E  1  =  e'J*04)  and  enclosed  by 
a  conformable  outer  boundary  as  shown  in  Figure  6.4.  Shown 
are  the  method  of  moments,  the  higher-order  boundary 
condition,  and  the  BGT  boundary  condition  results. 

a)  Real  part  of  the  near  field. 

b)  Imaginary  part  of  the  near  field. 
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Figure  6.6.  A  TM  incident  wave  on  a  6  X  p.e.c.  strip  enclosed  by  a  conformable 
outer  boundary  IX  away  to  minimize  the  number  of  mesh  points. 
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(b) 

Figure  6.8.  Near  field  on  the  outer  boundary  for  a  6?ip.e.c.  strip  illuminated  by 

a  TM  incident  wave  (E  *  =  e'J^)  and  enclosed  by  a  conformable 
outer  boundary  as  shown  in  Figure  6.7.  Shown  are  the  method  of 
moments,  the  higher-order  boundary  condition,  and  the  BGT 
boundary  condition  results. 

a)  Real  part  of  the  near  field. 

b)  Imaginary  part  of  the  near  field. 


Imaginary  part  of  near-field  Real  part  of  near-field 


(a) 


(b) 

Figure  6.9.  Same  as  Figure  6.8  except  in  this  graph  only  the  higher-order  and 
the  BGT  absorbing  boundary  condition  results  are  shown. 
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The  higher-order  boundary  condition  operator  depends  on  the  set  {ni,  n2,  n3}. 
This  set  is  chosen  in  a  way  that  includes  both  the  lower-  and  the  higher-order  harmonics. 
Let  y  denotes  the  ratio  up/u  (y  =  up/u  )  where  u  the  scattered  field  and  up  is  its  normal 
derivative.  The  imaginary  and  real  parts  of  y  are  related  to  the  propagating  and  decaying 
waves,  respectively.  The  exact  of  value  of  y  is  given  by  Hn’(kp)/Hn(kp).  The 
approximate  value  of  y  can  readily  by  computed  by  using  Equation  (3.5)  for  the  BGT 
boundary  condition  and  Equation  (6.12)  for  the  higher-order  boundary  condition.  Thus,  a 
proper  choice  of  the  harmonics  is  a  one  that  approximates  y  (Yapprox)  as  close  as  possible 
to  Yexact*  As  Figure  6.10  indicates  (kp  =  32),  the  BGT  boundary  condition  is  quite 
satisfactory  for  the  lower-order  harmonics  up  to  n=20,  and  then  begins  to  deviate  from  the 
exact  one  for  n>20.  However,  by  letting  the  set  {ni,  n2,  n3)  to  take  roughly  the  value  of 
{0,  kp/2,  kp},  where  k  is  the  free-space  wavenumber  and  p  is  the  distance  from  the  origin 
to  the  middle  of  the  triangular  edge  residing  on  the  outer  boundary,  the  higher-order 
absorbing  boundary  condition  yields  a  value  of  Yapprox  which,  on  the  average,  matches 
Yexact  on  a  wider  range  of  harmonics  than  does  the  one  given  by  the  BGT  boundary 
condition. 

As  the  numerical  results  indicate,  the  improvement  brought  about  by  the  higher- 
order  absorbing  boundary  condition  is  more  significant  for  the  wedge  than  for  the  strip. 
This  difference  can  be  explained  by  the  fact  that  the  strip  scatters  more  of  the  higher-order 
harmonics  than  does  the  wedge,  which,  in  turn,  may  require  more  than  three  sampling  of 
the  harmonics  in  order  to  achieve  a  significant  improvement 

The  term  ’’Error”  that  appeared  several  times  in  the  Tables  should  not  be  interpreted 
as  an  absolute  error  but  rather  as  a  difference  between  our  results  and  the  published  ones. 
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6.5  Conclusions 

Based  on  the  asymptotic  representation  of  the  general  solution  of  the  Helmholtz  and 
Laplace  equations,  we  derived  higher-order  absorbing  and  asymptotic  boundary  conditions 
that  combine  both  the  lower-  and  higher-order  terms.  The  higher-order  boundary 
conditions  were  then  generalized  and  made  valid  for  an  arbitrary  outer  boundary.  These 
boundary  conditions  were  implemented  in  the  finite  element  scheme  and  used  to  obtain  the 
solution  for  both  digital  circuit  and  scattering  applications.  The  numerical  results  showed 
that  the  higher-order  boundary  conditions  constantly  yielded  a  significant  improvement 
over  those  for  the  BGT  boundary  condition. 
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CHAPTER  7 

CONCLUSIONS  AND  FUTURE  WORK 

In  this  thesis,  the  subjects  of  the  absorbing  and  asymptotic  boundary  conditions  for 
the  finite  element  (FEM)  mesh  truncation  applied  to  scattering  and  digital  circuit  problems 
were  investigated.  Three  new  boundary  condition  concepts  were  introduced,  viz.,  the 
boundary  condition  for  arbitrary  outer  boundaries,  the  asymptotic  boundary  condition  for 
digital  circuit  applications,  and  the  higher-order  asymptotic  and  absorbing  boundary 
conditions. 

Using  a  rotation  and  a  translation  transformation  and  neglecting  the  mixed 
derivative  term  (un0.  the  most  commonly  used  absorbing  boundary  condition  operator  for 
electromagnetic  scattering  problems,  e.g.,  the  Bayliss,  Gunzburger,  and  Turkel  (BGT), 
was  generalized  and  made  applicable  to  an  arbitrary,  rather  than  circular,  outer  boundary 
for  the  purpose  of  minimizing  the  number  of  mesh  points.  The  generalized  boundary 
condition  operator  was  implemented  in  the  finite  element  method  and  used  to  study  the 
scattering  from  a  4X  strip,  a  2X  by  IX  wedge,  a  9X  strip,  and  an  8X  by  4X  wedge.  The 
numerical  investigation  indicated  that  while  the  finite  element  yielded  acceptable  results  for 
all  the  strips  and  wedges  considered,  the  results  for  the  wedges  agreed  more  favorably  than 
those  of  the  strips  with  the  results  obtainable  via  the  Method  of  Moments.  This  can  be 
explained  by  the  fact  that  the  scattered  waves  are  purely  outgoing  only  in  the  region  outside 
the  smallest  circle  that  entirely  encloses  the  scatterer.  For  a  wedge,  where  the  outer 
boundary  resembles  more  closely  a  circular  one,  there  are  more  points  satisfying  the  above 
criterion  than  there  are  for  a  strip  where  the  outer  boundary  is  elongated.  In  fact,  the 
presence  of  incoming  waves  is  more  pronounced  for  the  region  of  the  elongated  boundary 
enclosing  the  strip  which  is  close  to  the  origin..  The  issue  of  absorbing  boundary 
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condition  for  scattering  problems  was  examined  again  in  Chapter  6  where  a  higher-order 
boundary  condition  was  used  in  conjunction  with  the  finite  element  method  The  higher- 
order  absorbing  boundary  condition  was  derived  from  the  asymptotic  solution  of  the 
Helmholtz  equation.  Unlike  most  available  absorbing  boundary  conditions,  e.g.,  the 
Bayliss,  Gunzburger,  and  Turkel  (BGT),  'he  higher-order  absorbing  boundary  condition 
took  into  account  both  the  lower-  and  higher-order  harmonics.  Using  the  transformations 
of  Chapter  3,  this  boundary  condition  was  also  generalized  and  made  valid  for  an  arbitrary 
outer  boundary.  The  numerical  results  showed  that  the  higher-order  absorbing  boundary 
condition  constantly  yielded  a  significant  improvement  over  those  for  the  BGT  boundary 
condition  for  a  variety  of  scatterers.  In  the  course  of  this  study,  it  was  realized  that  whether 
one  uses  the  BGT  or  the  higher-order  boundary  conditions,  the  outer  boundary  needed  to 
be  placed  at  a  reasonable  distance  from  the  surface  of  the  scatterer  considered. 

The  corresponding  concept  to  the  absorbing  boundary  condition  for  digital  circuit 
applications  is  the  asymptotic  boundary  condition.  The  term  "asymptotic"  is  used  instead 
of  "absorbing"  because  the  digital  circuits  were  analyzed  in  the  quasi-TEM  regime  and, 
hence,  there  was  no  propagation  or  reflection  of  waves.  Using  a  similar  approach  to  the 
one  used  to  obtain  an  absorbing  boundary  condition  for  scattering  problems,  an  asymptotic 
boundary  condition  was  derived  for  open  region  digital  circuit  problems.  The  asymptotic 
boundary  condition  was  then  implemented  in  the  two-dimensional  FEM  scheme  to  model 
one-,  two-,  and  six-conductor  configurations.  The  numerical  results  indicated  that  the 
asymptotic  boundary  condition  consistently  yielded  more  accurate  results  than  those 
obtainable  with  a  perfectly  conducting  shield  placed  at  the  same  location.  Furthermore,  this 
asymptotic  boundary  condition  did  not  suffer  from  the  complications  associated  with  the 
infinite  elements.  However,  in  some  situations,  the  accuracy  obtained  with  the  asymptotic 
boundary  condition  was  not  adequate,  as  for  instance  in  the  off-diagonal  terms  of  the 
capacitance  matrix  in  the  six-conductor  configuration.  Those  inadequate  results  were 
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improved  in  Chapter  6  where  the  higher-order  asymptotic  boundary  condition  was  used. 
The  later  boundary  condition,  unlike  the  one  used  in  Chapter  3,  took  into  account  both  the 
lower-  and  higher-order  terms. 

Since  digital  circuits  do  not  consist  of  only  two-dimensional  structures  but  also 
contain  various  three-dimensional  transmission  line  discontinuities,  it  was  necessary  to 
derive  the  three-dimensional  version  of  the  asymptotic  boundary  condition.  The  general 
form  of  the  solution  to  Laplace's  equation  was  again  used  to  derive  the  asymptotic 
boundary  condition.  In  order  to  reduce  the  number  of  mesh  points  as  much  as  possible,  a 
box-shaped  outer  boundary  was  chosen  because  it  was  the  most  conformable  to  the 
structures  considered.  The  method  was  used  to  compute  the  capacitance  of  a  rectangular 
microstrip  patch,  and  the  results  were  found  to  be  in  good  agreement  with  data  published 
elsewhere.  Once  again,  the  asymptotic  boundary  condition  enabled  us  to  bring  the  outer 
boundary  much  closer  to  the  structure  than  would  have  been  possible  with  the  p.e.c. 
artificial  boundary.  Actually,  the  role  played  by  the  three-dimensional  asymptotic  boundary 
condition  was  more  crucial  than  the  one  played  by  the  two-dimensional  asymptotic 
boundary  condition  because  the  total  number  of  mesh  points  for  three-dimensional 
problems  is  usually  large.  Although  the  cost  of  FEM  calculation  for  the  three-dimensional 
problem  was  quite  high,  the  method  is  worth  using  because  it  handles  very  general 
structures  and,  hence,  helps  in  solving  practical  problems. 

The  absorbing  and  asymptotic  boundary  conditions  are  fairly  new  concepts  in  both 
scattering  and  digital  circuit  analyses.  The  work  presented  in  this  thesis  just  scratched  the 
surface  of  those  concepts  and  a  considerable  effort  still  needs  to  be  done  before  they 
become  standard  techniques  for  FEM  mesh  truncation.  A  possible  area  of  research  is 
suggested  by  the  fact  that  the  outer  boundary  had  to  be  placed  at  a  reasonable  distance  from 
the  object  in  order  to  achieve  an  acceptable  accuracy.  For  a  relatively  large  scatterer  such  as 
the  9X  strip  considered  earlier,  even  though  a  conformable  outer  boundary  was  used,  a 
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couple  of  thousand  mesh  points  were  required  in  order  to  obtain  satisfactory  results.  A 
considerable  savings  in  computer  time  and  storage  can  be  achieved  if  a  boundary  condition 
that  helps  to  bring  the  outer  boundary  closer  than  the  one  used  in  this  work  is  found.  But 
how  close  to  the  object  could  the  outer  boundary  be  placed?  The  answer  to  this  question  is 
still  an  unknown  issue  that  requires  a  great  deal  of  future  research. 

In  this  work,  only  two-dimensional  scattering  problems  were  considered. 
However,  most  of  the  practical  scattering  problems  are  three-dimensional  in  nature. 
Although  three-dimensional  vector  absorbing  boundary  condition  are  reported  in  the 
literature,  very  little  has  been  done  in  terms  of  implementing  them  in  the  FEM.  Thus,  there 
is  a  need  to  develop  computer  codes  that  incorporate  those  absorbing  boundary  condition 
operators  into  the  FEM  which  will  help  to  solve  the  most  general  scattering  problems,  i.e., 
the  vector  three-dimensional  problems. 

For  digital  circuit  problems,  the  asymptotic  boundary  conditions  derived  in  this 
work  were  valid  only  in  the  quasi-TEM  regime.  However,  with  the  increasing  speed  of 
digital  circuits,  it  is  often  inadequate  and  insufficient  to  rely  on  the  quasi-TEM 
approximations.  Thus,  there  is  an  urgent  need  to  derive  absorbing  boundary  conditions 
and  incorporate  them  into  the  FEM  for  the  full-wave  vector  analysis  of  two-  and  three- 
dimensional  high-speed  digital  circuit  problems. 
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